# Machine Learning-Based Estimation and Goodness-of-Fit for Large-Scale
Confirmatory Item Factor Analysis††thanks: This is a preprint of a paper
submitted to the journal Psychometrika.

Christopher J. Urban Correspondence to<EMAIL_ADDRESS>L.L. Thurstone
Psychometric Laboratory in the Department of Psychology and Neuroscience
University of North Carolina at Chapel Hill  Daniel J. Bauer  L.L. Thurstone
Psychometric Laboratory in the Department of Psychology and Neuroscience
University of North Carolina at Chapel Hill

###### Abstract

We investigate††This material is based upon work supported by the National
Science Foundation Graduate Research Fellowship under Grant No. DGE-1650116.
novel parameter estimation and goodness-of-fit (GOF) assessment methods for
large-scale confirmatory item factor analysis (IFA) with many respondents,
items, and latent factors. For parameter estimation, we extend Urban Bauer’s
(2021) deep learning algorithm for exploratory IFA to the confirmatory setting
by showing how to handle user-defined constraints on loadings and factor
correlations. For GOF assessment, we explore new simulation-based tests and
indices. In particular, we consider extensions of the classifier two-sample
test (C2ST), a method that tests whether a machine learning classifier can
distinguish between observed data and synthetic data sampled from a fitted IFA
model. The C2ST provides a flexible framework that integrates overall model
fit, piece-wise fit, and person fit. Proposed extensions include a C2ST-based
test of approximate fit in which the user specifies what percentage of
observed data can be distinguished from synthetic data as well as a C2ST-based
relative fit index that is similar in spirit to the relative fit indices used
in structural equation modeling. Via simulation studies, we first show that
the confirmatory extension of Urban Bauer’s (2021) algorithm produces more
accurate parameter estimates as the sample size increases and obtains
comparable estimates to a state-of-the-art confirmatory IFA estimation
procedure in less time. We next show that the C2ST-based test of approximate
fit controls the empirical type I error rate and detects when the number of
latent factors is misspecified. Finally, we empirically investigate how the
sampling distribution of the C2ST-based relative fit index depends on the
sample size.

Key words: Machine learning, deep learning, artificial neural network,
variational inference, item response theory, categorical factor analysis,
goodness-of-fit, fit indices

##  1 Introduction

Item factor analysis (IFA; Bock ., 1988) is an invaluable method for
investigating the latent structure underlying the discrete item response data
that arises in many social science applications. In particular, IFA allows
researchers to summarize a large number of item responses using a smaller
number of continuous latent factors, thereby reducing the dimensionality of
the data and potentially making the data easier to understand. Researchers
with specific hypotheses about the number of factors, the relations between
the item responses and the factors, and the factor correlations typically
encode their hypotheses as parameter constraints in a confirmatory IFA model
(e.g., Wirth Edwards, 2007). For example, personality assessments such as the
revised Minnesota Multiphasic Personality Inventory (Butcher ., 1989) and the
International Personality Item Pool NEO (Goldberg, 1999) include hundreds of
items organized into subscales wherein each subscale is designed to measure a
single personality factor. Modeling such designs using confirmatory IFA
entails estimating the relations between the items comprising a subscale and
their corresponding factor while constraining these items’ relations with all
other factors to zero (i.e., not estimating these relations). Estimating
confirmatory IFA model parameters and their standard errors permits inferences
about the properties of items as well as about the characteristics of the
population from which the observed sample was drawn. Subsequent goodness-of-
fit (GOF) analyses provide useful information about how well the fitted model
approximates the data generating model (e.g., Maydeu-Olivares, 20131).

Unfortunately, both parameter estimation and GOF assessment have long been
computationally challenging in the large-scale setting with many respondents,
items, and latent factors (e.g., Cai, 20102, 20101). Existing methods for
confirmatory IFA may therefore be sub-optimal for analyzing complex, high-
dimensional item response data arising from sources such as surveys,
standardized tests, online applications, and electronic data capture, many of
which are becoming increasingly available to social scientists (e.g., Pardos,
2017; Woo ., 2020). To clarify this issue, consider Bock Aitkin’s (1981)
marginal maximum likelihood (MML) estimator, which has many desirable
statistical properties and is typically the preferred estimator for
confirmatory IFA parameter estimation.111We note that other estimators such as
limited-information estimators (e.g., Jöreskog Moustaki, 2001; Muthén, 1978,
1984) and joint maximum likelihood (JML) estimators (X. Chen ., 2019) are more
computationally efficient than the MML estimator. However, these alternative
estimators have different statistical properties — for example, limited-
information estimators are not asymptotically efficient, while JML estimators
are only consistent when the sample size and the number of items
simultaneously tend to infinity — and are not considered further here due to
space constraints. The MML approach bases inference on the marginal likelihood
of the observed item responses, which is obtained by integrating out the
latent factors. Problematically, however, evaluating this integral is
computationally burdensome when the number of factors $P$ is even moderately
large (e.g., $P\geq$5$$). Researchers have devised numerous methods to avoid
this computational burden: adaptive Gaussian quadrature methods (Rabe-Hesketh
., 2005; Schilling Bock, 2005), Laplace approximation methods (e.g., Huber .,
2004), Monte Carlo expectation-maximization (EM) algorithms (e.g., Meng
Schilling, 1996; Song Lee, 2005), Markov Chain Monte Carlo methods (e.g.,
Béguin Glas, 2001; Edwards, 2010), and stochastic approximation methods (SA;
e.g., Cai, 20102, 20101; S. Zhang ., 2020). Of the above methods, SA
procedures such as the Metropolis-Hastings Robbins-Monro (MH-RM) algorithm
(Cai, 20102) and the stochastic EM (StEM) algorithm (S. Zhang ., 2020) are the
most computationally efficient; in recent years, MH-RM has been particularly
widely used in the social and behavioral sciences due to its flexibility and
computational efficiency. However, even these state-of-the-art SA procedures
are slow when the sample size $N$, number of items $J$, and number of factors
$P$ are all large (e.g., $N\geq$10\,000$$, $J\geq$100$$, and $P\geq$10$$).

Even after surmounting the computational difficulties associated with MML
estimation, researchers who wish to assess their confirmatory IFA model’s GOF
face yet another computational barrier. Let $K_{j}$ denote the number of
response categories for item $j$. To simplify the presentation, assume all
items have the same number of categories such that $K=K_{j}$ for
$j=1,\ldots,J$. GOF assessment for IFA models fitted via MML is typically
based on the underlying $K^{J}$-dimensional multinomial table on which the
model is defined. Full-information GOF statistics such as Pearson’s statistic
and the likelihood ratio statistic have inaccurate $p$-values when the number
of items and the number of response categories are even moderately large due
to the multinomial table’s sparseness (e.g., when $J\geq 6$ and $K\geq 5$;
Thissen Steinberg, 1997). Limited-information GOF statistics overcome the
sparsity problem by only utilizing marginals of the multinomial table, thereby
“concentrating” the information available for testing to obtain more accurate
p-values and higher power (e.g., Maydeu-Olivares Joe, 2005, 2006, 2014).
However, calculating limited-information GOF statistics requires high-
dimensional numerical integration and is typically computationally intensive.
Despite efforts to improve computational efficiency for specific IFA models
(e.g., Cai Hansen, 2013), calculating limited-information GOF statistics for
general confirmatory IFA models remains computationally intensive when the
number of items, response categories, and factors are all large (e.g., $J\geq
100$, $K\geq 5$, and $P\geq 10$).

Based on the preceding discussion, it is clear that more computationally
efficient MML estimation and GOF assessment methods are needed to apply
confirmatory IFA to very large-scale data. In this work, we investigate
machine learning methods that offer steps toward addressing some of the
difficulties mentioned above. The method we propose for parameter estimation
is based on Urban Bauer’s (2021) deep learning algorithm for exploratory IFA.
Their algorithm uses an importance-weighted amortized variational estimator
(I-WAVE) that combines variational inference and importance sampling to
construct an approximation to the MML estimator. By increasing the number of
importance-weighted samples drawn during fitting, the I-WAVE typically trades
computational efficiency for a better approximation. In the large-scale
exploratory setting (i.e., $P=10$, $J=100$, $$1000$\leq N\leq$10\,000$$), the
I-WAVE has empirically demonstrated comparable parameter estimation accuracy
and increased computational efficiency relative to the MML estimator
implemented via MH-RM (Urban Bauer, 2021).

Assessing GOF for large-scale IFA models fitted via I-WAVE is not
straightforward. In addition to being computationally inefficient in the
large-scale setting, the limited-information GOF statistics described above
were designed for models fitted via the MML estimator and have unknown
theoretical properties when applied to models fitted via approximate MML
(e.g., the I-WAVE). We aim to address both computational and theoretical
issues simultaneously by instead considering simulation-based GOF assessment
methods in which model-data fit is assessed by comparing the observed item
responses to synthetic item responses sampled from a fitted IFA model.
Previous work in simulation-based GOF assessment for confirmatory IFA has
mainly focused on posterior predictive model checking (PPMC) in the Bayesian
setting (e.g., Levy ., 2009; Sinharay ., 2006). In PPMC, synthetic data
simulated from the posterior predictive distribution are compared to the
observed data using discrepancy measures (i.e., measures that quantify how two
data sets differ). Unlike limited-information GOF statistics, PPMC is
computationally efficient, accounts for parameter estimation uncertainty, and
does not rely on asymptotic arguments. Although PPMC was originally developed
for Bayesian IFA models, recent work has shown that PPMC may be applied to
frequentist IFA models using a normal approximation to the posterior
predictive distribution (Lee ., 2016; Kuhfeld, 2019). Unfortunately, PPMC
assuming posterior normality (PPMC-N) is only well-motivated when models are
fitted via exact maximum likelihood, suggesting that PPMC-N may not be a well-
motivated GOF assessment method for I-WAVE.

The alternative GOF assessment methods we consider are based on a class of
machine learning methods called classifier two-sample tests (C2STs; Lopez-Paz
Oquab, 2017; Friedman, 2003; Kim ., 2021). Similar to $t$-tests and other
classical two-sample tests, C2STs aim to determine whether two samples are
drawn from the same distribution. To assess model-data fit in confirmatory
IFA, a C2ST begins by first sampling synthetic item responses from the fitted
IFA model.222Unlike PPMC-N, the C2ST is well-motivated when applied to an IFA
model fitted via any estimator given that synthetic data can be sampled from
the model. This holds for estimators that treat the latent factors as random
effects (e.g., the MML estimator or the I-WAVE) but not for estimators that
treat the latent factors as fixed effects (e.g., JML estimators). The
synthetic responses are combined with the observed responses to construct a
new data set, which is divided at random into two disjoint subsets called the
training set and the test set. Next, a machine learning classifier (e.g., a
logistic regression or a neural network classifier) is trained to distinguish
between the observed and synthetic training set response patterns. Finally,
the classifier’s test set accuracy (i.e., the proportion of test set response
patterns correctly classified as observed or synthetic) is treated as the test
statistic and is used to test the null hypothesis that the observed and
synthetic data are drawn from the same distribution. For a perfect-fitting
model, the accuracy obtained should not be significantly better than chance
because the synthetic data will mimic the observed data characteristics.
Depending on one’s choice for the classifier, C2STs are potentially very
computationally efficient in the large-scale setting. C2STs also provide a
variety of interpretable values to complement the use of $p$-values and may be
viewed as uniting overall model fit, piece-wise (i.e., item-level) fit, and
person fit under a single framework.

In exploring these issues, we make four primary contributions. First, we
extend the I-WAVE to the confirmatory setting by showing how to handle user-
defined constraints on the factor loadings and inter-factor correlations.
Second, we empirically investigate confirmatory I-WAVE’s finite sample
behavior and conduct comparisons with MH-RM. Third, in addition to a C2ST-
based test of perfect (exact) fit, we propose and explore a novel C2ST-based
test of approximate fit as well as a C2ST-based relative fit index that is
similar in spirit to the relative fit indices used in linear confirmatory
factor analysis and SEM (e.g., Bollen, 19892; Bentler Bonett, 1980; LR. Tucker
Lewis, 1973; Bentler, 1990). Compared to standard C2STs, these new methods are
potentially better suited to applications where the specified IFA model is
unlikely to exactly capture the data generating model. Fourth, we conduct
simulations to investigate the C2ST variants’ finite sample behavior in
several settings including when the IFA model is correctly and incorrectly
specified.

The remaining sections are organized as follows. Section 2 reviews variational
methods (e.g., I-WAVE) for fitting confirmatory IFA models with polytomous
responses. Novel C2ST-based GOF assessment methods are developed in Sect. 3.
Section 4 includes implementation details for the proposed parameter
estimation and GOF assessment methods as well as an empirical example and
simulation studies. Experimental results are summarized, then limitations and
extensions of the proposed methods are discussed in Sect. 5.

##  2 Variational Methods for Parameter Estimation

Variational inference (VI) is an approach to approximate maximum likelihood
estimation for latent variable (LV) models that is widely used in machine
learning (Blei ., 2017; C. Zhang ., 2019). VI has recently been applied for
IFA in a variety of settings (e.g., Curi ., 2019; Cho ., 2020; Y. Chen, Filho
., 2019; Natesan ., 2016; Wu ., 2020; Hui ., 2017; Urban Bauer, 2021). In this
section, we review variational methods for IFA. We focus in particular on the
VI-based method described by Urban Bauer (2021), which provides a theoretical
link between VI and MML estimation.

###  2.1 A Model for Confirmatory IFA

We first establish notation for the IFA models under consideration.
Specifically, we consider Samejima’s (1969) graded response model (GRM) for
polytomous item responses, although the methods discussed in this section
readily apply to other IFA models. Suppose that $N$ respondents have answered
$J$ items. Let $x_{i,j}\in\\{0,1,\ldots,K_{j}-1\\}$ denote the response for
respondent $i$ to item $j$ in $K_{j}$ ordinal categories. To simplify the
presentation, assume $K_{j}=K$ for $j=1,\ldots,J$. Note that when $K=2$, the
GRM reduces to the widely used multidimensional two-parameter logistic model
(McKinley Reckase, 1983).

Each respondent is represented by a $P\times 1$ latent vector $\mathbf{z}_{i}$
and each item is represented by a $(P+K-1)\times 1$ parameter vector
$\boldsymbol{\theta}_{j}=(\boldsymbol{\alpha}_{j}^{\top},\boldsymbol{\beta}_{j}^{\top})^{\top}$
where $\boldsymbol{\beta}_{j}$ is a $P\times 1$ vector of loadings and
$\boldsymbol{\alpha}_{j}=(\alpha_{j1},\ldots,\alpha_{j,K-1})^{\top}$ is a
$(K-1)\times 1$ vector of strictly ordered category intercepts. The GRM
defines a set of boundary response probabilities conditional on
$\boldsymbol{\theta}_{j}$ and $\mathbf{z}_{i}$:

$\Pr(x_{i,j}\geq
k\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i})=\frac{1}{1+\mathrm{exp}[\alpha_{j,k}+\boldsymbol{\beta}_{j}^{\top}\mathbf{z}_{i}]},\quad
k\in\\{1,\ldots,K-1\\},$ (1)

where $\Pr(x_{i,j}\geq 0\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i})=1$ and
$\Pr(x_{i,j}\geq K\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i})=0$. The
conditional probability of the response $x_{i,j}=k$ for $k\in\\{0,\ldots
K-1\\}$ is

$\pi_{i,j,k}=\Pr(x_{i,j}=k\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i})=\Pr(x_{i,j}\geq
k\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i})-\Pr(x_{i,j}\geq
k+1\mid\boldsymbol{\theta}_{j},\mathbf{z}_{i}).$ (2)

It follows from Equation (2) that the conditional distribution of $x_{i,j}$ is
multinomial with $K$ cells, trial size $1$, and cell probabilities
$\pi_{i,j,k}$:

$p_{\boldsymbol{\theta}_{j}}(x_{i,j}\mid\mathbf{z}_{i})=\prod_{k=0}^{K-1}\pi_{i,j,k}^{\mathbbm{1}(x_{i,j}=k)},$
(3)

where $\mathbbm{1}(\cdot)$ denotes the indicator function. Let
$\mathbf{x}_{i}=(x_{i,1},\ldots,x_{i,n})^{\top}$ be the $i^{\text{th}}$
respondent’s response pattern and let
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{1}^{\top},\ldots,\boldsymbol{\theta}_{J}^{\top})^{\top}$
be a vector collecting all item parameters. By the local independence
assumption, the conditional distribution of $\mathbf{x}_{i}$ is

$p_{\boldsymbol{\theta}}(\mathbf{x}_{i}\mid\mathbf{z}_{i})=\prod_{j=1}^{J}p_{\boldsymbol{\theta}_{j}}(x_{i,j}\mid\mathbf{z}_{i}).$
(4)

Assume that $\mathbf{z}_{i}$ is multivariate normally distributed with zero
mean vector and covariance matrix
$\boldsymbol{\Sigma}=(\sigma_{p,p^{\prime}})_{P\times P}$. Let
$\boldsymbol{\omega}=(\boldsymbol{\theta}^{\top},\text{vech}(\boldsymbol{\Sigma})^{\top})^{\top}$
be a vector collecting all unknown parameters where
$\text{vech}(\boldsymbol{\Sigma})$ stacks the $P(P+1)/2$ unique elements of
$\boldsymbol{\Sigma}$ into a vector. Under my assumptions about the
distribution of the factors, the marginal distribution of $\mathbf{x}_{i}$ is
given by

$p_{\boldsymbol{\omega}}(\mathbf{x}_{i})=\int\prod_{j=1}^{J}p_{\boldsymbol{\theta}_{j}}(x_{i,j}\mid\mathbf{z})\mathcal{N}(\mathbf{z}\mid\boldsymbol{\Sigma})d\mathbf{z},$
(5)

where $\mathcal{N}(\cdot\mid\boldsymbol{\Sigma})$ is a normal density
parameterized by $\boldsymbol{\Sigma}$ and the above integral is over
$\mathbb{R}^{P}$.

We set $\sigma_{p,p}=1$ for $p=1,\ldots,P$ to identify the scale of the
factors. In the confirmatory setting, users encode hypotheses about the
measurement structure by placing restrictions on the loadings (e.g., by fixing
$\beta_{j,p}$ to zero if item $j$ is not hypothesized to measure factor $p$).
Following Cai (20102), we consider the case of linear equality constraints so
that the loadings may be written as

$\boldsymbol{\beta}_{j}=\mathbf{b}_{j}+\mathbf{A}_{j}\boldsymbol{\beta}_{j}^{\prime},$
(6)

where $\boldsymbol{\beta}_{j}$ is the restricted loadings vector,
$\mathbf{b}_{j}$ is a $P\times 1$ vector of constants, $\mathbf{A}_{j}$ is a
$P\times P$ matrix of constants that implements the linear constraints, and
$\boldsymbol{\beta}_{j}^{\prime}$ is a vector of free parameters. Cai (20102)
provides examples of how $\mathbf{b}_{j}$ and $\mathbf{A}_{j}$ may be
specified to implement various restrictions, while Anderson Rubin (1956)
provide sufficient conditions enabling $\mathbf{b}_{j}$ and $\mathbf{A}_{j}$
to be specified such that the model is identified.

###  2.2 Variational Lower Bound

Let $\mathbf{X}$ be an $N\times J$ matrix whose $i^{\mathrm{th}}$ row is
$\mathbf{x}_{i}^{\top}$. The marginal log-likelihood of the observed data is

$\ell(\boldsymbol{\omega}\mid\mathbf{X})=\sum_{i=1}^{N}\log\bigg{[}\int\prod_{j=1}^{J}p_{\boldsymbol{\theta}_{j}}(x_{i,j}\mid\mathbf{z}_{i})\mathcal{N}(\mathbf{z}\mid\boldsymbol{\Sigma})d\mathbf{z})\bigg{]}.$
(7)

Maximizing $\ell(\boldsymbol{\omega}\mid\mathbf{X})$ by directly evaluating
the $N$ integrals in Equation (7) is computationally intensive when $P$ is
large. VI solves this issue by instead maximizing a computationally tractable
lower bound on $\ell(\boldsymbol{\omega}\mid\mathbf{X})$. To derive the lower
bound, we first re-write a single summand in Equation (7) as

$\log
p_{\boldsymbol{\omega}}(\mathbf{x}_{i})=\operatorname{\textit{D}_{KL}}\big{[}q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})\|p_{\boldsymbol{\omega}}(\mathbf{z}_{i}\mid\mathbf{x}_{i})\big{]}+\mathbb{E}_{q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})}\big{[}\log
p_{\boldsymbol{\omega}}(\mathbf{z}_{i},\mathbf{x}_{i})-\log
q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})\big{]},$ (8)

where $\operatorname{\textit{D}_{KL}}\big{[}\cdot\|\cdot\big{]}$ denotes the
Kullback-Leibler (KL) divergence333The KL divergence is defined as
$\operatorname{\textit{D}_{KL}}\big{[}\mathbb{Q}\|\mathbb{P}\big{]}=\mathbb{E}_{\mathbb{Q}}\big{[}\log\mathbb{Q}\big{]}-\mathbb{E}_{\mathbb{Q}}\big{[}\log\mathbb{P}\big{]}$
for distributions $\mathbb{P}$ and $\mathbb{Q}$. It satisfies
$\operatorname{\textit{D}_{KL}}\big{[}\mathbb{Q}\|\mathbb{P}\big{]}\geq 0$
with equality if and only if $\mathbb{P}=\mathbb{Q}$ almost everywhere. and
$q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})$ is an arbitrary density with
parameter vector $\boldsymbol{\psi}_{i}$ satisfying $\int
q_{\boldsymbol{\psi}_{i}}(\mathbf{z})d\mathbf{z}=1$. The first r.h.s. term in
Equation (8) is the KL divergence between
$q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})$ and the posterior distribution of
the latent factors. Since this term in nonnegative, the second r.h.s. term is
a lower bound on the marginal log-likelihood of a single observation:

$\log
p_{\boldsymbol{\omega}}(\mathbf{x}_{i})\geq\mathbb{E}_{q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})}\big{[}\log
p_{\boldsymbol{\omega}}(\mathbf{z}_{i},\mathbf{x}_{i})-\log
q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})\big{]}.$ (9)

This lower bound is called the evidence lower bound (ELBO) and is often re-
written as:

$\displaystyle\operatorname{ELBO}_{i}$
$\displaystyle=\mathbb{E}_{q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})}\big{[}\log
p_{\boldsymbol{\omega}}(\mathbf{z}_{i},\mathbf{x}_{i})-\log
q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})\big{]}$ (10)
$\displaystyle=\mathbb{E}_{q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})}\big{[}\log
p_{\boldsymbol{\omega}}(\mathbf{x}_{i}\mid\mathbf{z}_{i})\big{]}-\operatorname{\textit{D}_{KL}}\big{[}q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})\|\mathcal{N}(\mathbf{z}_{i}\mid\boldsymbol{\Sigma})\big{]}.$
(11)

The variational estimator of the IFA model parameters $\boldsymbol{\omega}$ is
obtained by maximizing the ELBO over all observations w.r.t. both
$\boldsymbol{\omega}$ and $\boldsymbol{\psi}_{i}$, which is equivalent to
minimizing the KL divergence between
$q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})$ and the posterior distribution of
the latent factors. Intuitively, obtaining the variational estimator drives
$q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})$ to approximate the true LV
posterior; we henceforth refer to $q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})$
as the approximate LV posterior. When the approximate LV posterior perfectly
approximates the true posterior such that the KL divergence term on the r.h.s.
of Equation (8) is zero, the ELBO equals the marginal log-likelihood and the
variational estimator is equivalent to the MML estimator.

Following previous work (e.g., Kingma Welling, 2014; Hui ., 2017; Urban Bauer,
2021), we set the approximate LV posterior to the isotropic normal density:

$q_{\boldsymbol{\psi}_{i}}(\mathbf{z}_{i})=\mathcal{N}(\mathbf{z}_{i}\mid\boldsymbol{\mu}_{i},\boldsymbol{\sigma}^{2}_{i}\mathbf{I}_{P}\big{)},$
(12)

where $\boldsymbol{\mu}_{i}$ is a $P\times 1$ vector of means,
$\boldsymbol{\sigma}^{2}_{i}$ is a $P\times 1$ vector of variances, and
$\mathbf{I}_{P}$ is a $P\times P$ identity matrix. Although the isotropic
normal approximate LV posterior is computationally tractable, it is rarely
flexible enough in practice to perfectly approximate the true posterior and
thereby minimize the KL divergence to zero. The importance sampling technique
described in Sect. 2.4 decreases the gap between the ELBO and the marginal
log-likelihood by implicitly increasing the flexibility of the approximate
posterior.

###  2.3 Amortized VI

Traditional VI fits a different parameter vector $\boldsymbol{\psi}_{i}$
(i.e., a different approximate LV posterior) for each observation. This
approach quickly becomes computationally infeasible for large sample sizes.
Additionally, it is not immediately clear how to apply models fitted using VI
to new observations (e.g., to compute the log-likelihood of the new
observations). Amortized variational inference (AVI) solves the above issues
by parameterizing the approximate posterior using a powerful function
approximator called an inference model. Since the inference model parameters
are shared across observations, performing AVI requires fitting a constant
number of parameters regardless of the sample size, whereas performing VI
requires fitting a number of parameters that grows linearly with the sample
size. Models fitted using AVI can also be straightforwardly applied to new
observations by passing the observations through the inference model.

The variational autoencoder (VAE; Kingma Welling, 2014; Rezende ., 2014) is an
algorithm for AVI that uses an artificial neural network (NN) inference model.
NNs are parametric machine learning models that map a set of predictor
variables through a sequence of transformations to predict a set of outcome
variables (for a brief overview, see Urban Bauer, 2021). Using an NN inference
model is typically considered justifiable because NNs can approximate any
Borel measurable function (e.g., Cybenko, 1989) and perform well in a variety
of real-world applications (LeCun ., 2015). We can specify a VAE for
confirmatory IFA by parameterizing the approximate LV posterior as follows:

$\displaystyle\begin{split}(\boldsymbol{\mu}_{i}^{\top},\log\boldsymbol{\sigma}_{i}^{\top})^{\top}&=\text{NN}_{\boldsymbol{\psi}}(\mathbf{x}_{i}),\\\
q_{\boldsymbol{\psi}}(\mathbf{z}_{i}\mid\mathbf{x}_{i})&=\mathcal{N}(\mathbf{z}_{i}\mid\boldsymbol{\mu}_{i},\boldsymbol{\sigma}^{2}_{i}\mathbf{I}_{P}),\end{split}$
(13)

where $\boldsymbol{\mu}_{i}$ is a $P\times 1$ predicted vector of means,
$\log\boldsymbol{\sigma}_{i}$ is a $P\times 1$ predicted vector of log-
standard deviations, and $\text{NN}_{\boldsymbol{\psi}}$ is a neural network
parameterized by $\boldsymbol{\psi}$. Instead of estimating a parameter vector
$\boldsymbol{\psi}_{i}$ for each observation, the NN parameters
$\boldsymbol{\psi}$ are now shared across observations. That is, instead of
maximizing Equation (11) over observations, we maximize

$\operatorname{ELBO}=\mathbb{E}_{q_{\boldsymbol{\psi}}(\mathbf{z}\mid\mathbf{x})}\big{[}\log
p_{\boldsymbol{\omega}}(\mathbf{x}\mid\mathbf{z})\big{]}-\operatorname{\textit{D}_{KL}}\big{[}q_{\boldsymbol{\psi}}(\mathbf{z}\mid\mathbf{x})\|\mathcal{N}(\mathbf{z}\mid\boldsymbol{\Sigma})\big{]}$
(14)

over observations. Note that we have dropped the case index $i$ since
$\boldsymbol{\psi}$ is shared across
$\\{\mathbf{x}_{i},\mathbf{z}_{i}\\}_{i=1}^{N}$.

The VAE relies on the assumption that the NN is capable of approximating the
mapping from the data space to the approximate LV posterior parameter space.
In theory, this assumption holds when the NN is sufficiently flexible (e.g.,
when the NN has one infinitely large hidden layer). In practice, however, the
NN has a limited capacity that may prevent the VAE from performing as well as
traditional VI. Fortunately, this performance difference can be reduced by
increasing the flexibility of the approximate LV posterior (Cremer ., 2018).

Technical details regarding model fitting are briefly reviewed here; see Urban
Bauer (2021) for greater detail. The fitting procedure of Urban Bauer (2021)
requires an unbiased estimator for the gradient of the ELBO w.r.t.
$\boldsymbol{\xi}=(\boldsymbol{\omega}^{\top},\boldsymbol{\psi}^{\top})^{\top}$.
This is obtained by first reparameterizing $\mathbf{z}$ as

$\displaystyle\begin{split}\boldsymbol{\epsilon}&\sim\mathcal{N}(\boldsymbol{\epsilon}),\\\
\mathbf{z}&=\boldsymbol{\mu}+\boldsymbol{\sigma}\odot\boldsymbol{\epsilon},\end{split}$
(15)

where $\boldsymbol{\epsilon}$ is a $P\times 1$ sample from a standard
multivariate normal density, $\boldsymbol{\mu}$ and $\boldsymbol{\sigma}$ are
the outputs of the NN inference model given in equations 13, and $\odot$
denotes element-wise multiplication. This reparameterization writes
$\mathbf{z}$ as a deterministic function of $\boldsymbol{\psi}$ and can be
viewed as “externalizing” the randomness in $\mathbf{z}$. An unbiased
estimator for the gradient of the ELBO w.r.t. $\boldsymbol{\xi}$ is now given
by

$\displaystyle\nabla_{\boldsymbol{\xi}}\operatorname{ELBO}$
$\displaystyle=\mathbb{E}_{\mathcal{N}(\boldsymbol{\epsilon})}\big{[}\nabla_{\boldsymbol{\xi}}\log
p_{\boldsymbol{\omega}}(\mathbf{z},\mathbf{x})-\nabla_{\boldsymbol{\xi}}\log
q_{\boldsymbol{\psi}}(\mathbf{z}\mid\mathbf{x})\big{]}$ (16)
$\displaystyle\approx\frac{1}{S}\sum_{s=1}^{S}\big{[}\nabla_{\boldsymbol{\xi}}\log
p_{\boldsymbol{\omega}}(\mathbf{z}_{s},\mathbf{x})-\nabla_{\boldsymbol{\xi}}\log
q_{\boldsymbol{\psi}}(\mathbf{z}_{s}\mid\mathbf{x})\big{]},$ (17)

where $\nabla_{\boldsymbol{\xi}}$ returns an $M\times 1$ vector of first-order
partial derivatives w.r.t. $\boldsymbol{\xi}$ and (17) is a Monte Carlo
approximation to the expectation in (16). Figure 1 illustrates how computation
proceeds in a VAE for confirmatory IFA: First, an item response $\mathbf{x}$
is passed through the inference model $\text{NN}_{\boldsymbol{\psi}}$ to
produce the approximate LV posterior parameters
$(\boldsymbol{\mu}^{\top},\log\boldsymbol{\sigma}^{\top})^{\top}$; next, a
latent vector $\mathbf{z}$ is sampled from
$\mathcal{N}(\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P}\big{)}$
and used to compute the response probabilities $\pi_{j,k}$ for $j=1,\ldots,6$
and $k=0,\ldots,2$; finally, these intermediate quantities are used to compute
both the expected conditional log-likelihood and the KL divergence terms in
the ELBO. The gradient estimator in Equation (17) can be efficiently computed
using an automatic differentiation procedure called backpropagation (e.g.,
Goodfellow ., 2016). Urban Bauer (2021) then apply an adaptive stochastic
gradient method called AMSGrad (Reddi ., 2018) to iteratively update
$\boldsymbol{\xi}$ until convergence.

$x_{1}$ $x_{2}$ $x_{3}$ $x_{4}$ $x_{5}$ $x_{6}$ $h_{1}$ $h_{2}$ $h_{3}$
$h_{4}$ $\mu_{1}$ $\mu_{2}$ $\log\sigma_{1}$ $\log\sigma_{2}$ Sample
$\mathbf{z}$ from
${\mathcal{N}\big{(}\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P}\big{)}}$
$z_{1}$ $z_{2}$ $\pi_{1,0}$ $\pi_{1,1}$ $\pi_{2,0}$ $\pi_{2,1}$ $\pi_{3,0}$
$\pi_{3,1}$ $\pi_{4,0}$ $\pi_{4,1}$ $\pi_{5,0}$ $\pi_{5,1}$ $\pi_{6,0}$
$\pi_{6,1}$
$\operatorname{\textit{D}_{KL}}\big{[}\mathcal{N}(\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P})\|\mathcal{N}(\mathbf{z}\mid\boldsymbol{\Sigma})\big{]}$$\underbrace{\phantom{\operatorname{\textit{D}_{KL}}\big{[}\mathcal{N}(\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P})\|\mathcal{N}(\mathbf{z}\mid\boldsymbol{\Sigma})\big{]}}}_{\text{Kullback-
Leibler
divergence}}$$\mathbb{E}_{\mathcal{N}(\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P})}\bigg{[}\sum\limits_{j=1}^{6}\sum\limits_{k=0}^{2}\mathbbm{1}(x_{j}=k)\log\pi_{j,k}\bigg{]}$$\underbrace{\phantom{\mathbb{E}_{\mathbf{z}\sim\mathcal{N}(\mathbf{z}\mid\boldsymbol{\mu},\boldsymbol{\sigma}^{2}\mathbf{I}_{P})}\bigg{[}\sum\limits_{j=1}^{6}\sum\limits_{k=0}^{2}\mathbbm{1}(x_{j}=k)\log\pi_{j,k}\bigg{]}}}_{\text{expected
conditional log-likelihood}}$ Input layer NN hidden layer Predicted
approximate LV posterior parameters LV layer Predicted response probabilities

Figure 1: Schematic diagram of a variational autoencoder for confirmatory item
factor analysis with $J=6$ items, $K=2$ categories per item, $P=2$ factors,
$S=1$ Monte Carlo sample from the approximate latent variable posterior, and
an inference model consisting of a neural network with a single hidden layer.
Each factor loads on three items. The reparameterization of $\mathbf{z}$ is
not illustrated for simplicity. LV = latent variable, NN = neural network.

###  2.4 Importance-Weighted VI

Importance-weighted VI (Burda ., 2016; Domke Sheldon, 2018) is a strategy for
obtaining a better approximation to the true log-likelihood by increasing the
flexibility of traditional VI. The importance-weighted amortized variational
estimator (I-WAVE) for the IFA model parameters $\boldsymbol{\omega}$ is
obtained by maximizing a new lower bound called the importance-weighted ELBO
(IW-ELBO):

$\displaystyle\log p_{\boldsymbol{\omega}}(\mathbf{x})$
$\displaystyle\geq\operatorname{IW-ELBO}$ (18)
$\displaystyle=\mathbb{E}_{\mathbf{z}_{1:R}}\bigg{[}\log\frac{1}{R}\sum_{r=1}^{R}w_{r}\bigg{]},$
(19)

where
$\mathbf{z}_{1:R}\sim\prod_{r=1}^{R}q_{\boldsymbol{\psi}}(\mathbf{z}_{r}\mid\mathbf{x})$,
$w_{r}=p_{\boldsymbol{\omega}}(\mathbf{z}_{r},\mathbf{x})/q_{\boldsymbol{\psi}}(\mathbf{z}_{r}\mid\mathbf{x})$
are unnormalized importance weights for the joint distribution of latent and
observed variables, and $R$ is the number of importance-weighted (IW) samples.
The IW-ELBO reduces to the ELBO when $R=1$ and converges monotonically to the
marginal log-likelihood as $R\rightarrow\infty$ under mild assumptions (Burda
., 2016). I-WAVE can also be viewed as implicitly defining a flexible
approximate LV posterior
$q_{\boldsymbol{\psi}}^{\text{IW}}(\mathbf{z}\mid\mathbf{x})$ that approaches
the true LV posterior pointwise as $R\rightarrow\infty$ (Cremer et al., 2017).
The above facts imply that I-WAVE and the MML estimator are equivalent when
the number of IW samples $R$ equals infinity, in which case I-WAVE inherits
the MML estimator’s statistical properties.

Following Urban Bauer (2021), we use Burda .’s (2016) estimator for the IW-
ELBO $\boldsymbol{\omega}$-gradient:

$\displaystyle\nabla_{\boldsymbol{\omega}}\mathbb{E}_{\mathbf{x}_{1:R}}\bigg{[}\log\frac{1}{R}\sum_{r=1}^{R}w_{r}\bigg{]}$
$\displaystyle=\mathbb{E}_{\boldsymbol{\epsilon}_{1:R}}\bigg{[}\sum_{r=1}^{R}\widetilde{w}_{r}\nabla_{\boldsymbol{\omega}}\log
w_{r}\bigg{]}$ (20)
$\displaystyle\approx\frac{1}{S}\sum_{s=1}^{S}\bigg{[}\sum_{r=1}^{R}\widetilde{w}_{r,s}\nabla_{\boldsymbol{\omega}}\log
w_{r,s}\bigg{]},$ (21)

as well as G. Tucker .’s (2019) “doubly reparameterized” estimator for the IW-
ELBO $\boldsymbol{\psi}$-gradient:

$\displaystyle\nabla_{\boldsymbol{\psi}}\mathbb{E}_{\mathbf{x}_{1:R}}\bigg{[}\log\frac{1}{R}\sum_{r=1}^{R}w_{r}\bigg{]}$
$\displaystyle=\mathbb{E}_{\boldsymbol{\epsilon}_{1:R}}\bigg{[}\sum_{r=1}^{R}\widetilde{w}_{r}^{2}\frac{\partial\log
w_{r}}{\partial\mathbf{z}_{r}}\frac{\partial\mathbf{z}_{r}}{\partial\boldsymbol{\psi}}\bigg{]}^{\top}$
(22)
$\displaystyle\approx\frac{1}{S}\sum_{s=1}^{S}\bigg{[}\sum_{r=1}^{R}\widetilde{w}_{r,s}^{2}\frac{\partial\log
w_{r,s}}{\partial\mathbf{z}_{r,s}}\frac{\partial\mathbf{z}_{r,s}}{\partial\boldsymbol{\psi}}\bigg{]}^{\top},$
(23)

where
$\boldsymbol{\epsilon}_{1:R}\sim\prod_{r=1}^{R}\mathcal{N}(\boldsymbol{\epsilon}_{r})$
and $\widetilde{w}_{r}=w_{r}/\sum_{r^{\prime}=1}^{R}w_{r^{\prime}}$ are
normalized importance weights. Both estimators are unbiased and have
increasing signal-to-noise ratios as $R\rightarrow\infty$. Further, both can
be successfully approximated using a single Monte Carlo sample (e.g., Burda .,
2016; G. Tucker ., 2019), so we set $S=1$ for all experiments in this work.

###  2.5 Handling User-Defined Constraints

User-defined constraints on the factor loadings are straightforward to
implement for I-WAVE. In particular, the gradient of the IW-ELBO w.r.t. the
unconstrained loadings vector $\boldsymbol{\beta}_{j}^{\prime}$ can be
obtained using the chain rule:

$\nabla_{\boldsymbol{\beta}_{j}^{\prime}}\operatorname{IW-
ELBO}=\bigg{(}\frac{\partial\operatorname{IW-
ELBO}}{\partial\boldsymbol{\beta}_{j}}\frac{\partial\boldsymbol{\beta}_{j}}{\partial\boldsymbol{\beta}_{j}^{\prime}}\bigg{)}^{\top}=\mathbf{A}^{\top}_{j}\nabla_{\boldsymbol{\beta}_{j}}\operatorname{IW-
ELBO}.$ (24)

Equation (24) implies that one can first compute the gradient of the IW-ELBO
w.r.t. $\boldsymbol{\beta}_{j}$, then obtain the gradient w.r.t.
$\boldsymbol{\beta}_{j}^{\prime}$ via pre-multiplication by the transposed
constraint matrix $\mathbf{A}^{\top}_{j}$.

Users also often wish to impose constraints on the factor correlation matrix
$\boldsymbol{\Sigma}$. Let $\boldsymbol{\Sigma}=\mathbf{L}\mathbf{L}^{\top}$
where $\mathbf{L}$ is a $P\times P$ lower triangular matrix. We estimate
$\mathbf{L}$ using a hyperspherical parameterization (Pinheiro Bates, 1996;
Rapisarda ., 2007), which enables unconstrained estimation of a variety of
structured correlation matrices and has similar computational efficiency to
estimating $\mathbf{L}$ directly. This parameterization is given by:

$l_{p,p^{\prime}}=\begin{cases}\cos\vartheta_{p,1},&\text{if}\;p^{\prime}=1\\\
\cos\vartheta_{p,p^{\prime}+1}\prod_{p^{\prime\prime}=1}^{p^{\prime}}\sin\vartheta_{p,p^{\prime\prime}},&\text{if}\;1<p^{\prime}<p\\\
\prod_{p^{\prime\prime}=1}^{p}\sin\vartheta_{p,p^{\prime\prime}},&\text{if}\;p^{\prime}=p,\end{cases}$
(25)

for $p=1,\ldots,P$ where $l_{p,p^{\prime}}$ are elements of $\mathbf{L}$ and
$\vartheta_{p,p^{\prime}}\in(0,\pi]$ are angles measured in radians which are
elements of a $P\times P$ lower triangular matrix $\boldsymbol{\Theta}$.
Constraints on the angles giving rise to various correlation structures are
discussed by Tsay Pourahmadi (2017) as well as by Ghosh . (2020). We note that
it is also feasible to estimate $\boldsymbol{\Sigma}$ using a proximal
stochastic gradient method (e.g., Yun ., 2020).

##  3 Classifier Two-Sample Tests for Goodness-of-Fit Assessment

###  3.1 Exact C2STs

We now discuss the application of classifier two-sample tests (C2STs;
Friedman, 2003; Kim ., 2021; Lopez-Paz Oquab, 2017), a class of simulation-
based GOF assessment methods that have recently been developed in machine
learning, to assessing exact GOF for confirmatory IFA models. Let
$\hat{\boldsymbol{\omega}}$ denote parameter estimates obtained for some
confirmatory IFA model. Let $\mathbf{x}_{i}\sim\mathbb{P}$ denote the
$i^{\text{th}}$ observed response pattern and let $\mathbf{y}_{j}\sim
p_{\hat{\boldsymbol{\omega}}}(\mathbf{y}_{j})=\hat{\mathbb{P}}$ denote the the
$j^{\text{th}}$ synthetic response pattern drawn from the fitted model where
$\mathbf{x}_{i},\mathbf{y}_{j}\in\mathcal{X}$ for $i=1,\ldots,N_{1}$ and
$j=1,\ldots,N_{2}$. To simplify the presentation, we assume $N_{1}=N_{2}=N$.

C2STs aim to test whether the observed and synthetic response patterns are
drawn from the same distribution — that is, C2STs aim to test
$H_{0}:\mathbb{P}=\hat{\mathbb{P}}$ against
$H_{1}:\mathbb{P}\neq\hat{\mathbb{P}}$. A C2ST is conducted by training a
machine learning classifier to distinguish between the observed and synthetic
response patterns. Intuitively, when $\mathbb{P}=\hat{\mathbb{P}}$, the
classifier’s test set accuracy should be close to $1/2$ (i.e., chance), since
samples from $\mathbb{P}$ and $\hat{\mathbb{P}}$ are indistinguishable. When
$\mathbb{P}\neq\hat{\mathbb{P}}$, the classifier should be able to capitalize
on the distributional differences to obtain a test set accuracy higher than
$1/2$. More formally, a C2ST can be described as the following five-step
procedure (Lopez-Paz Oquab, 2017):

  1. (1)

Construct a data set
$D=\\{(\mathbf{x}_{i},1)\\}_{i=1}^{N}\cup\\{(\mathbf{y}_{i},0)\\}_{i=1}^{N}=\\{(\mathbf{u}_{i},l_{i})\\}_{i=1}^{2N}$.

  2. (2)

Shuffle $D$ at random and split it into disjoint sets
$D=D_{\operatorname{train}}\cup D_{\operatorname{test}}$ where
$N_{\operatorname{train}}=|D_{\operatorname{train}}|$ and
$N_{\operatorname{test}}=|D_{\operatorname{test}}|$.

  3. (3)

Fit a machine learning classifier $\hat{f}:\mathcal{X}\rightarrow[0,1]$ using
$D_{\operatorname{train}}$.

  4. (4)

Obtain the test set classification accuracy:

$\widehat{\operatorname{\textit{acc}}}=\frac{1}{N_{\operatorname{test}}}\sum_{i\in
I_{\operatorname{test}}}\mathbbm{1}\left(\mathbbm{1}\left(\hat{f}(\mathbf{u}_{i})>\frac{1}{2}\right)=l_{i}\right)=\frac{1}{N_{\operatorname{test}}}\sum_{i\in
I_{\operatorname{test}}}\widehat{\operatorname{\textit{acc}}}_{i},$ (26)

where $I_{\operatorname{test}}=\\{i:(\mathbf{u}_{i},l_{i})\in
D_{\operatorname{test}}\\}$.

  5. (5)

Compute a $p$-value based on the asymptotically normal null distribution of
$\operatorname{\textit{acc}}$:

$\hat{p}=\Pr(\widehat{\operatorname{\textit{acc}}}^{\prime}\geq\widehat{\operatorname{\textit{acc}}}\mid
H_{0})=1-\Phi\left(\frac{\widehat{\operatorname{\textit{acc}}}-1/2}{\sqrt{\frac{1}{4N_{\operatorname{test}}}}}\right),$
(27)

where $\Phi(\cdot)$ denotes the standard normal cumulative distribution
function.

To derive Equation (27), notice that
$\operatorname{\textit{acc}}_{i}\sim\operatorname{Bernoulli}(\operatorname{\textit{acc}}_{i}\mid
p_{i})$ where $p_{i}=1/2$ is the probability of correctly classifying some
$\mathbf{u}_{i}$ in the test set when $H_{0}$ is true. In this setting,
Equation (27) follows from the fact that

$N_{\operatorname{test}}\operatorname{\textit{acc}}\sim\operatorname{Binomial}\left(N_{\operatorname{test}}\operatorname{\textit{acc}}\,\middle|\,N_{\operatorname{test}},\frac{1}{2}\right)\approx\mathcal{N}\left(N_{\operatorname{test}}\operatorname{\textit{acc}}\,\middle|\,\frac{N_{\operatorname{test}}}{2},\frac{N_{\operatorname{test}}}{4}\right)$
(28)

as $N_{\operatorname{test}}\rightarrow\infty$.

C2STs provide a variety of interpretable numbers that complement the use of
$p$-values:

  1. (a)

Taking $\hat{f}(\mathbf{u}_{i})$ as an estimate of the conditional probability
$\Pr(l_{i}=1\mid\mathbf{u}_{i})$ for $i\in I_{\operatorname{test}}$, we can
determine which item response patterns were labeled correctly or incorrectly
as well as how confident $\hat{f}$ was in each decision. This approach
provides a way to evaluate which observed response patterns are discrepant
from the fitted IFA model.

  2. (b)

We can interpret the fitted classifier $\hat{f}$ to determine which items were
most useful for distinguishing between real and synthetic distributions (e.g.,
by interpreting a logistic regression classifier’s coefficients).

  3. (c)

We can interpret the test statistic $\widehat{\operatorname{\textit{acc}}}$ as
the percentage of item response patterns that were correctly distinguished
between the real and synthetic distributions.

The values described in (a), (b), and (c) correspond to measures of person
fit, piece-wise fit, and overall model fit, respectively. C2STs may therefore
be viewed as uniting these different kinds of fit measures under a single
framework.

C2STs are closely related to classical two-sample tests including the $t$-test
for the difference in means of two samples (Student, 1908), the Wilcoxon-Mann-
Whitney test for the difference in rank means of two samples (Wilcoxon, 1945;
Mann Whitney, 1947), and the Kolmogorov-Smirnov test for the difference in
empirical cumulative distributions of two samples (Kolmogorov, 1933; Smirnov,
1939), although these classical tests are less statistically efficient in the
multivariate setting. C2STs are also related to modern multivariate kernel
two-sample tests (K2STs) including the maximum mean discrepancy test (Gretton
., 2012) and the mean embedding test (Jitkrittum ., 2016; Chwialkowski .,
2015), both of which test for differences in the empirical kernel mean
embeddings of two samples. However, depending on one’s choice for $f$, C2STs
typically improve upon K2STs by requiring less manual engineering of the input
data and by providing more interpretable results.

###  3.2 Approximate C2STs

The C2STs described above are exact in the sense that they test the null
hypothesis that the real distribution $\mathbb{P}$ and the synthetic
distribution $\hat{\mathbb{P}}$ are exactly equal. In general, however, it is
unlikely that any specified IFA model will exactly capture the data generating
mechanism such that $\mathbb{P}=\hat{\mathbb{P}}$ (e.g., Cudeck Henly, 1991;
Maccallum Tucker, 1991). We therefore propose a more realistic approximate
C2ST for which we assume $\mathbb{P}\neq\hat{\mathbb{P}}$ and we test
$H_{0}:\operatorname{\textit{acc}}=1/2+\delta$ against
$H_{1}:\operatorname{\textit{acc}}>1/2+\delta$ where $\delta\in(0,1/2)$ is a
pre-specified value representing the degree of model error viewed as tolerable
by the user. The approximate C2ST is not a test of exact GOF because it does
not test whether the hypothesized IFA model exactly captures the data
generating mechanism. Instead, the approximate C2ST is a test of approximate
GOF wherein the user asserts that an IFA model that fits the data “well
enough” should be capable of synthesizing item response patterns that can only
be distinguished from real item response patterns around $100\delta\%$ of the
time.

We now derive the null distribution of $\operatorname{\textit{acc}}$ for the
approximate C2ST. Since $\operatorname{\textit{acc}}=1/2+\delta>1/2$ under
$H_{0}$, Equation (26) implies that
$\operatorname{\textit{acc}}_{i}=1/2+\delta_{i}\geq 1/2$ for $i\in
I_{\operatorname{test}}$ where $\delta_{i}\in[0,1/2]$ and
${\delta=N_{\operatorname{test}}^{-1}\sum_{i\in
I_{\operatorname{test}}}\delta_{i}}$. In this setting, the
$\operatorname{\textit{acc}}_{i}$ are independent but not identically
distributed Bernoulli random variables with success probabilities
$p_{i}=1/2+\delta_{i}$. $N_{\operatorname{test}}\operatorname{\textit{acc}}$
therefore follows a Poisson binomial distribution, which we follow Ehm (1991)
in approximating as

$N_{\operatorname{test}}\operatorname{\textit{acc}}\sim\operatorname{Binomial}\left(N_{\operatorname{test}}\operatorname{\textit{acc}}\mid
N_{\operatorname{test}},\bar{p}\right)\approx\mathcal{N}\big{(}N_{\operatorname{test}}\operatorname{\textit{acc}}\mid
N_{\operatorname{test}}\bar{p},N_{\operatorname{test}}\bar{p}(1-\bar{p})\big{)}$
(29)

as $N_{\operatorname{test}}\rightarrow\infty$ where
$\bar{p}=N_{\operatorname{test}}^{-1}\sum_{i\in
I_{\operatorname{test}}}p_{i}=1/2+N_{\operatorname{test}}^{-1}\sum_{i\in
I_{\operatorname{test}}}\delta_{i}=1/2+\delta$. It follows that

$\operatorname{\textit{acc}}\sim\mathcal{N}\left(\operatorname{\textit{acc}}\,\middle|\,\frac{1}{2}+\delta,\frac{\frac{1}{4}-\delta^{2}}{N_{\operatorname{test}}}\right).$
(30)

To derive the alternative distribution of $\operatorname{\textit{acc}}$,
notice that under $H_{1}$ we have
$\operatorname{\textit{acc}}=1/2+\delta+\varepsilon>1/2$ where the effect size
$\varepsilon\in(0,1/2-\delta)$ is the magnitude of the difference between
$\operatorname{\textit{acc}}$ under $H_{0}$ and $\operatorname{\textit{acc}}$
under $H_{1}$. Equation (26) implies that
$\operatorname{\textit{acc}}_{i}=1/2+\delta_{i}+\varepsilon_{i}$ for $i\in
I_{\operatorname{test}}$ where $\delta_{i}\in[0,1/2]$,
$\varepsilon_{i}\in[0,1/2-\delta_{i}]$,
${\delta=N_{\operatorname{test}}^{-1}\sum_{i\in
I_{\operatorname{test}}}\delta_{i}}$, and
${\varepsilon=N_{\operatorname{test}}^{-1}\sum_{i\in
I_{\operatorname{test}}}\varepsilon_{i}}$. Then by a similar argument to the
one given in the previous paragraph, we can obtain

$\operatorname{\textit{acc}}\sim\mathcal{N}\left(\operatorname{\textit{acc}}\,\middle|\,\frac{1}{2}+\delta+\varepsilon,\frac{\frac{1}{4}-\delta^{2}-2\delta\varepsilon-\varepsilon^{2}}{N_{\operatorname{test}}}\right)$
(31)

as $N_{\operatorname{test}}\rightarrow\infty$.

We now analyze the approximate C2ST’s power (i.e., the probability of
correctly rejecting $H_{0}$ when $H_{0}$ is false) by proving the following
theorem.

######  Theorem 3.1.

Let $\alpha\in[0,1]$ be the user-defined significance level (i.e., the
probability of incorrectly rejecting $H_{0}$ when $H_{0}$ is true). Suppose
the null and alternative distributions of $\operatorname{\textit{acc}}$ are
given by equations 30 and 31, respectively. Then the power of the approximate
C2ST is approximately given by

$\operatorname{power}(\alpha,N_{\operatorname{test}},\delta,\varepsilon)=\Phi\left(\frac{\varepsilon\sqrt{N_{\operatorname{test}}}-\sqrt{\frac{1}{4}-\delta^{2}}\Phi^{-1}(1-\alpha)}{\sqrt{\frac{1}{4}-\delta^{2}-2\delta\varepsilon-\varepsilon^{2}}}\right).$

Proof of Theorem 3.1 is given in Appendix A and follows the approach of Lopez-
Paz Oquab (2017, Theorem 1).

######  Remark 1.

The (approximate) power of the exact C2ST was derived by Lopez-Paz Oquab
(2017, Theorem 1) and can be alternately be derived by setting $\delta=0$ in
the power formula in Theorem 3.1.

Theorem 3.1 provides an important insight: maximizing the approximate C2ST’s
power involves trading off between maximizing the test set size
$N_{\operatorname{test}}$ and choosing a classifier $f$ that maximizes the
effect size $\varepsilon$ (equivalently, that maximizes the test set accuracy
$\operatorname{\textit{acc}}=1/2+\delta+\varepsilon$). To illustrate this
trade-off, notice that while relatively inflexible classifiers (e.g., logistic
regression) often obtain low test set accuracy because they cannot capture
complex nonlinear relationships (i.e., they decrease $\varepsilon$), these
classifiers only require small training data sets (i.e., they increase
$N_{\operatorname{test}}$). Conversely, while flexible classifiers (e.g., NNs)
may obtain high test set accuracy by capturing complex nonlinear relationships
(i.e., they increase $\varepsilon$), these classifiers require large training
data sets (i.e., they decrease $N_{\operatorname{test}}$). Lopez-Paz Oquab
(2017) note that this trade-off is analogous to the well known bias-variance
trade-off in statistics and machine learning (e.g., Yarkoni Westfall, 2017).
In this work, we empirically investigate this tradeoff by using two different
classifiers (i.e., a $k$-nearest neighbors classifier and an NN classifier)
that may require different training set sizes to obtain optimal test set
accuracy. We also follow Jitkrittum . (2016) and Lopez-Paz Oquab (2017) in
setting $N_{\operatorname{train}}=N_{\operatorname{test}}=N$, which often
achieves high power in practice and would achieve maximum power if
$\mathbb{P}$ and $\hat{\mathbb{P}}$ differed only in means.

###  3.3 C2ST-Based Relative Fit Index

The C2STs described above measure how well the proposed IFA model reproduces
the observed data (i.e., the model’s absolute GOF). An alternative approach
that is often used in SEM is to calculate a relative fit index (RFI) that
measures the proportional improvement in fit obtained by moving to the
proposed model from a more restrictive baseline model (e.g., Bentler Bonett,
1980; Bentler, 1990; Bollen, 19891; LR. Tucker Lewis, 1973). A typical
baseline model posits that the observed variables are mutually independent
(i.e., there are no common latent factors underlying the data; Bentler Bonett,
1980) and serves as a contrasting point of reference to a model that perfectly
reproduces the observed data (i.e., a saturated model; e.g., Steiger, 1980;
Bentler, 1995).

Consider a zero-factor baseline model where $P=0$ and
$\boldsymbol{\theta}_{j}=\boldsymbol{\alpha}_{j}$ for $j=1,\ldots,J$. In this
case, each respondent’s response probability $\pi_{i,j,k}$ evaluated at the
maximum likelihood estimate of $\boldsymbol{\theta}$ can be shown to be the
observed proportion of respondents choosing response category $k$ for item
$j$, which we write as
$\hat{\pi}_{j,k}=N^{-1}\sum_{i=1}^{N}\mathbbm{1}(x_{i,j}=k)$. We can therefore
conduct a C2ST for this baseline model using $N$ synthetic samples drawn from
a multinomial distribution with $K$ cells, trial size $1$, and cell
probabilities $\hat{\pi}_{j,k}$. Let
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$ and
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$ denote the test
set accuracies for the proposed model and for the baseline model,
respectively, and let $\hat{f}_{\operatorname{prop}}$ and
$\hat{f}_{\operatorname{base}}$ denote the corresponding fitted classifiers.
We propose the following C2ST-based RFI:

$\operatorname{C2ST-
RFI}=1-\frac{M_{\operatorname{prop}}}{M_{\operatorname{base}}}\cdot\frac{\Delta_{\operatorname{prop}}}{\Delta_{\operatorname{base}}},$
(32)

where
$\Delta_{\operatorname{prop}}=\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}-1/2$,
$\Delta_{\operatorname{base}}=\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}-1/2$,
$M_{\operatorname{prop}}$ is the number of fitted parameters in the proposed
model, and $M_{\operatorname{base}}$ is the number of fitted parameters in the
baseline model. The ratio
$\Delta_{\operatorname{prop}}/\Delta_{\operatorname{base}}$ is a measure of
the proportional change in misfit obtained by moving to the proposed model
from the baseline model, while the ratio
$M_{\operatorname{prop}}/M_{\operatorname{base}}$ is a penalty that increases
with the number of fitted parameters in the proposed model (i.e., it rewards
parsimony in the proposed model).

We motivate our definition of $\operatorname{C2ST-RFI}$ by considering the fit
index’s behavior in the typical setting where the baseline model obtains less-
than-perfect fit (i.e., when
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}>1/2$). In the
common scenario that the proposed model fits the same or better than the
baseline model (i.e., when
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$ varies between
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$ and $1/2$),
$\operatorname{C2ST-RFI}$ varies between
$1-M_{\operatorname{prop}}/M_{\operatorname{base}}$ and one with values closer
to one indicating better fit. $\operatorname{C2ST-RFI}$ is much larger than
one when $\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$ is much
smaller than $1/2$, which may occur when proposed model fits the data well but
$\hat{f}_{\operatorname{prop}}$ has overfitted the training data. In the
infrequent event that the proposed model fits worse than the baseline model
(i.e., when
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}>\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$),
$\operatorname{C2ST-RFI}$ is smaller than
$1-M_{\operatorname{prop}}/M_{\operatorname{base}}$. $\operatorname{C2ST-RFI}$
therefore behaves somewhat analogously to other nonnormed fit indices such as
the Tucker-Lewis Index (LR. Tucker Lewis, 1973) that distinguish between less-
than-perfect model fit, perfect model fit, and overfitting by being smaller
than one in the first case, close to one in the second case, and much larger
than one in the third case.

We next consider the effect of sample size on $\operatorname{C2ST-RFI}$.
Bollen (19891) notes that GOF indices may be influenced by sample size either
(a) when $N$ enters the calculation of the index or (b) when the mean of the
sampling distribution of the index is related to $N$. For (a), notice that
although both $\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$
and $\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$ include
factors of $N_{\operatorname{test}}^{-1}$, these factors cancel in the ratio
$\Delta_{\operatorname{prop}}/\Delta_{\operatorname{base}}$, indicating that
(a) does not hold. For (b), notice that when $N$ is very small, both
$\hat{f}_{\operatorname{prop}}$ and $\hat{f}_{\operatorname{base}}$ may fail
to capture any relationships in the training data. In this case, both
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$ and
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$ will be close to
$1/2$ and $\operatorname{C2ST-RFI}$ will be close to
$1-M_{\operatorname{prop}}/M_{\operatorname{base}}$. As $N$ increases,
$\hat{f}_{\operatorname{prop}}$ and $\hat{f}_{\operatorname{base}}$ should
capture any relationships in the training data and
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$ and
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}$ should come
close to their population values, suggesting that (b) holds. We provide
empirical evidence that the mean of the sampling distribution of
$\operatorname{C2ST-RFI}$ depends on $N$ in Sect. 4. Specifically, we observe
that C2ST-RFI gives a less optimistic assessment of fit as $N$ increases. This
phenomenon is analogous to the approximate C2ST’s power increasing with $N$
and is likely of little concern for moderate- to large-scale applications.

We also consider the probability limit of the proposed fit index. Assume that
$\operatorname{plim}_{N\rightarrow\infty}(\operatorname{\textit{acc}}_{\operatorname{prop}})=1/2+\delta$
and
$\operatorname{plim}_{N\rightarrow\infty}(\operatorname{\textit{acc}}_{\operatorname{base}})=1/2+\delta+\varepsilon$
where $\delta\in[0,1/2]$ and $\varepsilon\in[0,1/2-\delta]$. This assumption
is mild and states that the accuracies obtained by
$\hat{f}_{\operatorname{prop}}$ and $\hat{f}_{\operatorname{base}}$ should
come close to specific constants as $N$ grows large, with
$\hat{f}_{\operatorname{prop}}$ obtaining the same or better accuracy than
$\hat{f}_{\operatorname{base}}$.444A similar assumption that
$\operatorname{plim}_{N,\;J\rightarrow\infty}(\operatorname{\textit{acc}})=1/2+\delta$
for some $\delta>0$ is used by Kim . (2021) to prove that C2STs are consistent
(i.e., have power approaching one) in the high-dimensional setting where both
$N$ and $J$ tend to infinity. Then
$\operatorname{plim}_{N\rightarrow\infty}(\operatorname{C2ST-
RFI})=1-[M_{\operatorname{prop}}/M_{\operatorname{base}}]\cdot[\delta/(\delta+\varepsilon)]=c$
where $c$ is a constant which equals one when the proposed model is correct
(i.e., when $\delta=0)$ and is smaller than one otherwise.

Finally, we comment on how to interpret numerical values of C2ST-RFIs. In
practice, cutoff values are used to interpret fit indices and thereby evaluate
model fit. The development of adequate “rules of thumb” cutoff criteria that
cover a wide range of realistic data and model conditions requires extensive
empirical study (e.g., Hu Bentler, 1999) and is beyond the scope of this work.
Based on limited initial experimentation, however, we consider a provisory
cutoff of $\operatorname{C2ST-RFI}>0.9$ to indicate good fit. We empirically
investigate the performance of this cutoff in Sect. 4.

##  4 Implementation and Experiments

###  4.1 Implementation

I-WAVE is programmed using the machine learning library PyTorch (Version 1.1.6
Paszke ., 2017). Although PyTorch supports GPU computing to accelerate
fitting, we follow Urban Bauer (2021) in opting for CPU computing to enable
fairer comparisons with other methods and to assess performance using hardware
that is more typically available to social scientists. Experiments are
conducted on a computer with a 2.8 GHz Intel Core i7 CPU and 16 GB of RAM.
Code to reproduce all experiments is readily available at
https://github.com/cjurban/MachineLearningConfirmatoryIFA.

Optimization and hyperparameter tuning for I-WAVE in the exploratory setting
are discussed in detail by Urban Bauer (2021, Section 5). We follow their
implementation almost exactly, although we adapt their factor loadings
initialization procedure to the confirmatory setting as follows:

$\beta_{j,p}\sim\operatorname{Uniform}\left(\beta_{j,p}\,\middle|\,-\sqrt{\frac{6}{M_{p}+P}},\sqrt{\frac{6}{M_{p}+P}}\right),$
(33)

for $j=1,\ldots,J$ and $p=1,\ldots,P$ where $M_{p}$ is the number of nonzero
loadings parameters for factor $p$. This strategy is based on Xavier
initialization (Glorot Bengio, 2010) and aims stabilize fitting while taking
into account user-defined zero constraints on the factor loadings. An
additional detail that must be addressed in the confirmatory setting is
initializing the hyperspherical parameterization $\boldsymbol{\Theta}$ of the
factor correlation matrix $\boldsymbol{\Sigma}$. To ensure that
$\boldsymbol{\Sigma}$ is well-conditioned at the start of fitting, we choose
$\boldsymbol{\Theta}$ such that $\boldsymbol{\Sigma}=\mathbf{I}_{P}$ by
setting $\vartheta_{p,p^{\prime}}=\pi/2$ for $p=1,\ldots,P$,
$p^{\prime}=1,\ldots,p-1$.

Following Lopez-Paz Oquab (2017), we investigate two classifiers for C2ST
analyses: a $k$-nearest neighbors (KNN) classifier (Fix Hodges, 1951) and an
NN classifier.555We initially also investigated a regularized logistic
regression classifier but found that it performed very poorly relative to KNNs
and NNs. Both classifiers are implemented via the Python package scikit-learn
(Version 0.23.1; Pedregosa ., 2011) and have different hyperparameter
settings, which we discuss in turn. For each test set observation
$\mathbf{u}_{i}\in\mathcal{X}$, let $I_{i}$ collect the indices of the $k$
training set observations that are closest to $\mathbf{u}_{i}$ under some
distance metric $d:\mathcal{X}\times\mathcal{X}\rightarrow[0,\infty)$. KNNs
set the predicted class $\hat{l}_{i}$ for $\mathbf{u}_{i}$ to the most
frequent observed class for the $k$ closest training set observations, that
is, $\hat{l}_{i}=\mathbbm{1}(k^{-1}\sum_{i^{\prime}\in
I_{i}}l_{i^{\prime}}>1/2)$. In this work, we use the Hamming distance metric
$d(\mathbf{u}_{i},\mathbf{u}_{i^{\prime}})=J^{-1}\sum_{j=1}^{J}\mathbbm{1}(u_{i,j}\neq
u_{i^{\prime},j})$, which returns the proportion of unequal responses between
$\mathbf{u}_{i}$ and $\mathbf{u}_{i^{\prime}}$. We also set
$k=\lfloor\sqrt{N_{\operatorname{test}}}\rfloor$, which often performs
adequately in practice (Duda ., 2001, Chapter 4).

Although KNNs are flexible and often require very little tuning, they
typically obtain worse predictive accuracy as $J$ increases and are
computationally intensive when the training set size
$N_{\operatorname{train}}$ is large (Hastie ., 2009, Chapter 13). NNs overcome
these issues in that their predictive accuracy suffers relatively less as $J$
increases and their computational efficiency does not depend on
$N_{\operatorname{train}}$ when they are fitted via stochastic gradient
methods (Bottou ., 2018), although they often require more tuning than KNNs to
obtain optimal predictive accuracy. We set most NN hyperparameters to the
scikit-learn package defaults, which performed well in all experiments. To
mitigate overfitting in small data sets, we (a) stop fitting if the default
convergence criterion is not obtained within $\lfloor$10\,000$\cdot
200/N_{\operatorname{test}}\rfloor$ stochastic gradient steps and (b) apply
weight decay (i.e., an $L_{2}$ norm penalty on the NN regression weights)
where the optimal weight decay hyperparameter is selected from
$\\{10^{n}:n\in\\{-1,-1/2,0,1/2,1\\}\\}$ using a held-out validation set
consisting of $25\%$ of the training set observations (for details, see Hastie
., 2009, Chapters 7 and 11).

Since neither KNNs nor NNs are directly interpretable, we use a permutation
importance (PI) method to interpret which items are most useful for
distinguishing between real and synthetic distributions (Breiman, 2001; Fisher
., 2019). PI can be used to interpret any fitted classifier $\hat{f}$, even
when $\hat{f}$ itself is not directly interpretable. Let $\mathbf{U}$ be an
$N\times J$ matrix whose $i^{\mathrm{th}}$ row is $\mathbf{u}^{\top}_{i}$ for
$i\in I_{\operatorname{test}}$. For each item $j=1,\ldots,J$ and each
repetition $t=1,\ldots,T$, shuffle the $j^{\mathrm{th}}$ column of
$\mathbf{U}$ to generate a corrupted matrix $\widetilde{\mathbf{U}}_{j,t}$
whose $i^{\mathrm{th}}$ row is $\widetilde{\mathbf{u}}^{\top}_{i,j,t}$. The PI
for item $j$ is defined as the average decrease in test set accuracy for
$\hat{f}$ when the responses to item $j$ are shuffled:

$\operatorname{\textit{importance}}_{j}=\widehat{\operatorname{\textit{acc}}}-\frac{1}{T}\sum_{t=1}^{T}\frac{1}{N_{\operatorname{test}}}\sum_{i\in
I_{\operatorname{test}}}\mathbbm{1}\bigg{(}\mathbbm{1}\bigg{(}\hat{f}(\widetilde{\mathbf{u}}_{i,j,t})>\frac{1}{2}\bigg{)}=l_{i}\bigg{)}=\widehat{\operatorname{\textit{acc}}}-\frac{1}{T}\sum_{t=1}^{T}\widetilde{\operatorname{\textit{acc}}}_{j,t}.$
(34)

Since the shuffling procedure breaks the relationship between item $j$ and the
class label, $\operatorname{\textit{importance}}_{j}$ measures how much
$\hat{f}$ depends on item $j$. PI may be viewed as a piece-wise fit assessment
method wherein each $\operatorname{\textit{importance}}_{j}$ measures how well
the proposed IFA model fits item $j$, with larger
$\operatorname{\textit{importance}}_{j}$ indicating worse fit. We set the
number of repetitions to $T=5$ for all experiments.

###  4.2 Empirical Example

We demonstrate the proposed methods’ computational efficiency in the large-
scale setting and obtain data generating parameters for simulation studies by
analyzing $1\,015\,342$ responses to the $50$-item International Personality
Item Pool five-factor model (IPIP-FFM; Goldberg, 1999), the same data
considered by Urban Bauer (2021) for exploratory IFA using I-WAVE. The IPIP-
FFM is designed to measure the Big Five personality factors of openness,
conscientiousness, extraversion, agreeableness, and emotional stability,
making this data well suited to a more confirmatory approach. Each factor is
measured by $10$ five-category items anchored by “Disagree” (1), “Neutral”
(3), and “Agree” (5). The data were downloaded from the Open-Source
Psychometrics Project (https://openpsychometrics.org/) and pre-processed
following Urban Bauer (2021, Section 6.1), resulting in an analytic sample
containing $N=$515\,708$$ responses. Reverse worded items were recoded so that
the highest numerical response category indicated a high level of the
corresponding factor.

####  4.2.1 A Five-Factor Model

We used I-WAVE to fit a five-factor confirmatory IFA model with correlated
factors to the IPIP-FFM data. Optimization and inference model hyperparameters
were set to the values used by Urban Bauer (2021, Section 6.1). We set the
number of IW samples to $R=5$ based on Urban Bauer’s (2021) finding that this
value performs well in practice. We fitted the data set $10$ times to
investigate whether parameter estimation and GOF assessment results were
stable across random starts.

Factor loadings and correlations from the fitted model that attained the
highest IW-ELBO across random starts — henceforth called the reference model —
are reported in Tables 1 and 2, respectively. All loadings were positive,
which fit with the confirmatory design of the measurement scale. Factor
correlations aligned with the typical finding that emotional stability is
negatively correlated with the other factors. The mean fitting time was $198$
seconds ($SD=56$ seconds), which is quite fast given the large sample size.
Relative to the reference model, mean loadings root-mean-square error (RMSE)
was $0.031$ ($SD=0.017$), mean intercepts RMSE was $0.046$ ($SD=0.018$), and
mean factor correlation RMSE was $0.027$ ($SD=0.010$), suggesting that fitting
was fairly stable.666We treat the model attaining the highest IW-ELBO as
ground truth in line with how an optimal solution is often selected from
multiple random starts for mixture models (e.g., Biernacki ., 2003). Were we
to alternatively treat the mean parameter estimates across random starts as
ground truth, the corresponding variability estimates would be smaller than
those reported here, although this approach would downplay the impact of
possible local maxima on parameter estimate stability.

Table 1: Factor Loadings for IPIP-FFM Data Set Factor | Loading | Items
---|---|---
Five-Factor Model
Extraversion | $2.12$ | $2.18$ | $2.05$ | $2.59$ | $2.57$ | $1.55$ | $2.38$ | $1.64$ | $1.73$ | $2.27$ | 1–10
Emotional Stability | $2.09$ | $1.47$ | $1.69$ | $1.22$ | $1.27$ | $2.35$ | $2.35$ | $2.66$ | $1.92$ | $1.94$ | 11–20
Agreeableness | $1.54$ | $1.72$ | $.74$ | $2.80$ | $2.16$ | $1.45$ | $2.17$ | $1.59$ | $2.24$ | $1.14$ | 21–30
Conscientiousness | $1.69$ | $1.53$ | $.83$ | $1.80$ | $1.73$ | $1.95$ | $1.39$ | $1.36$ | $1.56$ | $.97$ | 31–40
Openness | $1.32$ | $1.59$ | $1.50$ | $1.43$ | $1.71$ | $1.65$ | $1.14$ | $1.17$ | $.85$ | $2.30$ | 41–50
Seven-Factor Model
Extraversion | $2.00$ | $2.15$ | $2.06$ | $2.57$ | $2.62$ | $1.54$ | $2.38$ | $1.56$ | $1.74$ | $2.33$ | 1–10
Emotional Stability | $2.40$ | $1.59$ | $1.90$ | $1.23$ | $1.38$ | $2.39$ | $2.87$ | $3.27$ | $1.94$ | $1.83$ | 11–20
Agreeableness | $1.53$ | $1.74$ | $.70$ | $2.85$ | $2.19$ | $1.43$ | $2.18$ | $1.56$ | $2.21$ | $1.12$ | 21–30
Conscientiousness | $1.78$ | $1.49$ | $.87$ | $1.77$ | $1.75$ | $1.92$ | $1.34$ | $1.35$ | $1.49$ | $1.03$ | 31–40
Openness | $1.69$ | $1.48$ | $1.56$ | $1.38$ | $1.80$ | $1.77$ | $1.03$ | $1.43$ | $.83$ | $2.54$ | 41–50
Doublet 1 |  |  |  |  |  |  | $2.31$ | $2.31$ |  |  | 11–20
Doublet 2 | $2.38$ |  |  |  |  |  |  | $2.38$ |  |  | 41–50
Table 2: Factor Correlations for IPIP-FFM Data Set | Factor
---|---
Factor | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$
Five-Factor Model
1\. Extraversion | $1.00$ |  |  |  |  |  |
2\. Emotional Stability | $-.23$ | $1.00$ |  |  |  |  |
3\. Agreeableness | $.37$ | $.01$ | $1.00$ |  |  |  |
4\. Conscientiousness | $.10$ | $-.28$ | $.08$ | $1.00$ |  |  |
5\. Openness | $.21$ | $-.10$ | $.18$ | $.01$ | $1.00$ |  |
Seven-Factor Model
1\. Extraversion | $1.00$ |  |  |  |  |  |
2\. Emotional Stability | $-.26$ | $1.00$ |  |  |  |  |
3\. Agreeableness | $.38$ | $.02$ | $1.00$ |  |  |  |
4\. Conscientiousness | $.09$ | $-.27$ | $.11$ | $1.00$ |  |  |
5\. Openness | $.26$ | $-.15$ | $.20$ | $.06$ | $1.00$ |  |
6\. Doublet 1 | $.00$ | $.00$ | $.00$ | $.00$ | $.00$ | $1.00$ |
7\. Doublet 2 | $.00$ | $.00$ | $.00$ | $.00$ | $.00$ | $.00$ | $1.00$

We assessed overall fit for the proposed five-factor model (FFM) across random
starts using KNN- and NN-based exact and approximate C2STs. The large sample
size required minor modifications to how the classifiers were implemented.
Specifically, we fitted the KNNs using a random subsample of $2.5\%$ of the
training set observations to ensure that fitting was not computationally
intensive, and we did not apply weight decay to the NNs given that large
sample sizes often mitigate overfitting. For approximate C2STs, we set
$\delta=0.025$ to test $H_{0}:\operatorname{\textit{acc}}=0.525$ against
$H_{1}:\operatorname{\textit{acc}}>0.525$, which amounts to testing whether
the proposed FFM is capable of synthesizing item response patterns that can
only be distinguished from real item response patterns $52.5\%$ of the time
(i.e., at slightly better than chance). C2ST results are presented in Table 3.
The fitted KNNs obtained much lower test set accuracies than the fitted NNs
($M=0.67$ for KNNs vs. $M=0.83$ for NNs), echoing the typical finding that NNs
outperform KNNs when $J$ is large. All KNN- and NN-based C2STs rejected
$H_{0}$ at signficance level $\alpha=0.05$ ($\hat{p}<0.0001$), suggesting that
the FFM did not fit the data “well enough”. Importantly, all tests were fast:
accounting for both sampling synthetic data and fitting the classifier, KNN-
based C2STs took around $20$ seconds and NN-based C2STs took around $15$
seconds.

Table 3: Classifier Two-Sample Test Results for IPIP-FFM Data Set

|  | Baseline Model |  | Five-Factor Model |  | Seven-Factor Model
---|---|---|---|---|---|---
|  | $M$ | $SD$ |  | $M$ | $SD$ |  | $M$ | $SD$
Sampling Time |  | $6.8$ |  $<$$.1$ |  | $8.7$ | $.4$ |  | $8.7$ | $.7$
KNN |  |  |  |  |  |  |  |  |
Fitting Time |  | $11.0$ | $.7$ |  | $11.8$ | $2.2$ |  | $11.6$ | $.7$
Test Set Accuracy |  | $.65$ | $.01$ |  | $.67$ | $.01$ |  | $.66$ | $.01$
C2ST-RFI |  |  |  |  | $-.49$ | $.22$ |  | $-.46$ | $.21$
NN |  |  |  |  |  |  |  |  |
Fitting Time |  | $6.5$ |  $<$$.1$ |  | $6.6$ | $.3$ |  | $6.6$ |  $<$$.1$
Test Set Accuracy |  | $.96$ |  $<$$.01$ |  | $.83$ |  $<$$.01$ |  | $.81$ |  $<$$.01$
C2ST-RFI |  |  |  |  | $.05$ | $.02$ |  | $.10$ | $.01$

  * •

Note. “Sampling Time” refers to time required to sample synthetic data from
the model. All times are given in seconds.

We next assessed piece-wise fit by computing PIs for fitted classifiers
accross random starts, which are displayed in Figure 2(a). The median NN
importances are clearly much larger than the median KNN importances, likely
due to the NNs having obtained higher test set accuracies. Further, the NNs
rely heavily on five specific items — items $17$, $18$, and $20$, which
measure emotional stability, and items $41$ and $48$, which measure openness —
whereas the KNNs rely fairly evenly on all the items. This finding suggests
that the NNs captured relationships between these five items and the class
label that the KNNs did not, likely contributing to the NNs’ superior
performance. The NN results in particular suggest that the FFM’s overall lack
of fit is most heavily influenced by its piece-wise lack of fit to items $17$,
$18$, $20$, $41$, and $48$.

(a) Five-factor model.

(b) Seven-factor model.

Figure 2: Permutation importances for each IPIP-FFM item across $10$ random
starts. Items are grouped by their corresponding factor. Bar heights indicate
medians, while error bars indicate $25\%$ and $75\%$ quantiles. Hatched bars
indicate the five poorest fitting items for the five-factor model. EXT =
extraversion, EST = emotional stability, AGR = agreeableness, CON =
conscientiousness, OPN = openness.

Although the above results suggest that the FFM did not fit the data well in
an absolute sense, similar findings might be expected for nearly any a priori
model. It would therefore be useful to assess the FFM’s value in a relative
sense by contrasting it with a baseline model. To this end, we computed
${\operatorname{C2ST-RFI}}\mathrm{s}$ for the fitted classifiers across random
starts. As shown in Table 3, KNNs obtained much lower baseline model test set
accuracies than NNs ($M=0.65$ for KNNs vs. $M=0.96$ for NNs). The KNN results
suggest that the proposed FFM and the baseline zero-factor model had
comparable fit, while the NN results suggest that the FFM fit the data better
than the baseline model. This finding was reflected in the
${\operatorname{C2ST-RFI}}\mathrm{s}$, which were much lower for KNNs than for
NNs ($M=-0.49$ for KNNs vs. $M=0.05$ for NNs). No C2ST-RFIs exceeded the
provisory cutoff of $0.9$ suggested as an indicator of good fit. The
additional RFI computations were also fast, taking at most around $18$ seconds
or less.

####  4.2.2 A Seven-Factor Model

We investigated the five poorest fitting items flagged by PIs with the goal of
improving model fit. As shown in Table 4, wordings for item pair $17$ and $18$
as well as for item pair $41$ and $48$ are similar. To account for possible
local dependence between these similarly worded pairs, we modeled each pair
using an additional orthogonal “doublet” factor that was only measured by its
corresponding two items and whose loadings were constrained to be equal to
ensure model identification. The resulting seven-factor model (SFM) was fitted
with $10$ random starts using the same hyperparameters as the FFM.

Table 4: Wordings for Five Poorest Fitting IPIP-FFM Items Item | Wording
---|---
17 | I change my mood a lot.
18 | I have frequent mood swings.
20 | I often feel blue.
41 | I have a rich vocabulary.
48 | I use difficult words.

The SFM loadings and factor correlation estimates are given in Tables 1 and 2,
respectively, and are largely similar to the FFM estimates for the non-doublet
factors. Fitting remained fast ($M=188$ seconds, $SD=44$ seconds) and fairly
stable (relative to the seven-factor reference model, loadings RMSE $M=0.038$,
$SD=0.016$; intercepts RMSE $M=0.051$, $SD=0.021$; and factor correlation RMSE
$M=0.032$, $SD=0.010$).

C2ST results for the SFM in Table 3 suggest that overall fit improved
marginally relative to the FFM. In particular, mean test set accuracy for both
KNNs and NNs was slightly lower for the SFM, suggesting that the SFM fit the
data slightly better. However, all exact and approximate C2STs again rejected
$H_{0}$ at $\alpha=0.05$ ($\hat{p}<0.0001$), suggesting that the SFM also
failed to fit the data “well enough”. The SFM’s fit relative to the zero-
factor baseline also improved only marginally: mean KNN- and NN-based C2ST-
RFIs were slighly higher for the SFM than for the FFM, although these values
remained far from the provisory cutoff of $0.9$. Computation for all tests and
fit indices remained fast.

Although the SFM appeared to improve overall fit only marginally relative to
the FFM, piece-wise fit appeared to improve more substantially. This
improvement is evident in the SFM’s PIs, which are shown in Figure 2(b). In
particular, NN-based PIs for the flagged item pairs are drastically lower for
the SFM than for the FFM. Improvements in piece-wise fit for these item pairs
were also observed for KNNs, although these improvements were relatively small
since KNN-based PIs were much smaller than NN-based PIs for both the FFM and
the SFM.

###  4.3 Simulation Studies

####  4.3.1 Evaluating I-WAVE

We investigate confirmatory I-WAVE’s parameter recovery and computational
efficiency as the number of IW samples $R$ increases and the log-likelihood
approximation improves. The data generating model has $P=5$ factors and $J=50$
$5$-category items. Generating parameters are rounded estimates from the five-
factor reference model in Section 4.2.1. We consider $R=1$, $5$, and $25$ as
well as $N=500$, $2500$, $12\,500$, and $62\,500$, resulting in $12$ total
simulation settings for each $R$ and $N$ combination. We conduct $100$
replications at each setting. All analyses reused the optimization and
inference model hyperparameters from Section 4.2.

Parameter recovery was assessed by computing the bias for each parameter as
the mean deviation of the estimated parameter from the data generating
parameter across replications, that is,
$\operatorname{bias}(\hat{\xi},\xi)=100^{-1}\sum_{a=1}^{100}[\hat{\xi}^{(a)}-\xi]$
where $\hat{\xi}^{(a)}$ is the estimated parameter at replication $a$ and
$\xi$ is the data generating parameter. We also computed mean squared error
(MSE) for each parameter as
$\operatorname{MSE}(\hat{\xi},\xi)=100^{-1}\sum_{a=1}^{100}[\hat{\xi}^{(a)}-\xi]^{2}$.
Boxplots of parameter biases and MSEs for each simulation setting are
displayed in Figures 3 and 4, respectively, with separate plots for factor
loadings, factor correlations, and intercepts. All estimates become more
accurate as the sample size increases. Loadings and factor correlations are
highly biased for $R=1$ but appear nearly unbiased for $R\geq 5$. MSE tends to
decrease with increasing $R$ for each $N$ setting. We note that several
replications initially converged to poor local maxima when $N=500$ and $R\geq
5$ (specifically, $7$ replications for $R=5$ and $10$ replications for
$R=25$). When $N\geq$2500$$ and $R\geq 5$, poor local maxima occurred in one
or fewer replications at each $(R,N)$ setting. Problematic replications were
handled by refitting with a new random seed.

Figure 3: Parameter bias for the importance-weighted amortized variational
estimator (I-WAVE). Three settings for the number of importance-weighted (IW)
samples are compared. Figure 4: Parameter mean squared error (MSE) for I-WAVE.

Line plots of fitting times for each simulation setting are displayed in
Figure 5. Median fitting time decreases from around $116$ seconds to around
$96$ seconds as $R$ increases from $1$ to $5$, then increases to around $165$
seconds as $R$ increases from $5$ to $25$. For fixed $R$, median fitting time
tends to decrease slightly as $N$ increases. These results demonstrate that
I-WAVE is computationally efficient even when the sample size is very large.

Figure 5: Fitting times for I-WAVE. Markers indicate medians, while error bars
indicate $25\%$ and $75\%$ quantiles.

####  4.3.2 Comparing I-WAVE to MH-RM

In this study, we compare I-WAVE to the MML estimator implemented via MH-RM in
a setting where the number of factors is large. We use the MH-RM
implementation from the R package mirt, which has core functions written in
both R and C++ (Version 1.32.1 Chalmers, 2012). The data generating model has
$P=10$ factors measured by $J=100$ $5$-category items. Generating parameters
are again rounded estimates from the five-factor reference model in Section
4.2.1 with the parameters for items $51$–$100$ set equal to the parameters for
items $1$–$50$. The factor correlation matrix is a $10\times 10$ block
diagonal matrix with rounded FFM estimates on the main-diagonals and zeros on
the off-diagonals. We conduct $100$ replications for each $N=$1000$$, $2500$,
$5000$, and $10\,000$. For I-WAVE, we set $R=5$ and set the learning rate
optimization hyperparameter to $0.0025$ to ensure the IW-ELBO does not diverge
due to the large number of factors (for details, see Urban Bauer, 2021).
I-WAVE converged to poor local maxima once each when $N=$1000$$ and $2500$,
both of which were handled by refitting with a new random seed. We set MH-RM
hyperparameters to the mirt package defaults, which performed well across $N$
settings.

Simulation results are presented in Figures 6 and 7. Both methods obtain
comparable estimates in all $N$ settings and have MSE decreasing toward zero
with increasing $N$. I-WAVE may obtain slightly better estimates than MH-RM
for $N=$1000$$, and vice versa for $N=$10\,000$$, although the differences
appear to be somewhat negligible. Parameter biases for both methods were also
comparable and are not shown. I-WAVE is faster than MH-RM in all settings: the
median fitting time for I-WAVE stays around three minutes and decreases
slightly as $N$ increases, whereas the median fitting time for MH-RM is around
$4.5$ minutes when $N=$1000$$ and increases to around $18$ minutes when
$N=$10\,000$$. We note that the MH-RM implementation in the commercially
available flexMIRT software (Cai, 2017) has core functions written in C++ and
is likely faster than mirt, although even this implementation would become
slower as $N$ increases. Additionally, I-WAVE could achieve a potentially
large speedup using a GPU, whereas a similar speedup is not currently
available for MH-RM.

Figure 6: MSE for I-WAVE and the marginal maximum likelihood estimator. MH-RM
= Metropolis-Hastings Robbins-Monro. Figure 7: Fitting times for I-WAVE and
the marginal maximum likelihood estimator.

####  4.3.3 Evaluating Approximate C2STs in a Non-IFA Setting

It is challenging to verify that approximate C2STs have accurate empirical
type I error rates and power in the confirmatory IFA setting because precisely
controlling the effect size $\varepsilon$ is not straightforward. We instead
verify these properties using real-valued data drawn from tractable
distributions whereby $\varepsilon$ is precisely controlled. Descriptions and
results of these experiments are provided in Appendix B. To briefly summarize,
we find that empirical type I error rates are well controlled and that
empirical power stays close to the values predicted by the formula in Theorem
3.1.

####  4.3.4 Evaluating C2STs and C2ST-RFIs in a Confirmatory IFA Setting

We now investigate the proposed GOF assessment methods’ performance in
settings where the fitted IFA model is correctly or incorrectly specified. We
consider two data generating models for $J=50$ five-category items: (1) the
same FFM considered in Section 4.3.1 as well as (2) an SFM with five
correlated factors each measured by $10$ items, one orthogonal doublet factor
measured by items $17$ and $18$, and another orthogonal doublet factor
measured by items $41$ and $48$. Generating parameters for (2) are rounded
estimates from the seven-factor reference model in Section 4.2.2. We simulate
$100$ data sets from (1) and (2) for each $N=$750$$, $1250$, $2500$, $5000$,
and $10\,000$. Using I-WAVE with the same hyperparameters as in Section 4.2,
we fit two models to each data set: (a) an FFM with the same specification as
(1) and (b) an SFM with the same specification as (2). This leads to four
different settings for each combination of data generating model and fitted
model. Viable GOF assessment methods should indicate near perfect fit for
settings (1a) and (2b), which have correctly specified fitted models, as well
as for setting (1b), which has an overspecified fitted model. Setting (2a), on
the other hand, has an underspecified fitted model and should demonstrate
poorer fit as $N$ increases. Although it might be expected that
$\operatorname{\textit{acc}}=0.5$ in settings with correctly specified and
overspecified models, it is also feasible that IFA model parameter estimate
uncertainty leads to $\operatorname{\textit{acc}}$ slightly greater than
$0.5$. To investigate both possibilities, we assess GOF for all simulation
settings using KNN- and NN-based exact C2STs where $\delta=0$ as well as
approximate C2STs where $\delta=0.025$, which respectively correspond to
$H_{0}:\operatorname{\textit{acc}}=0.5$ and
$H_{0}:\operatorname{\textit{acc}}=0.525$.

Rejection rates at signficance level $\alpha=0.05$ as well as test set
classification accuracies for the overspecified SFM in setting (1b) are shown
in Figure 8(a). Results for the correctly specified models in settings (1a)
and (2b) were nearly identical and are not shown. Rejection rates remain well
below the nominal level and test set classification accuracies come close to
$0.5$ as $N_{\operatorname{test}}$ increases, suggesting that the SFM fits the
simulated data nearly perfectly. The near-zero rejection rates likely occurred
because there was almost no signal in the training data, leading classifiers
to overfit to noise and perform worse than random chance for small to moderate
$N$. A similar phenomenon was observed in Appendix B, where it appeared to
have a negligible impact on empirical power.

(a) An overspecified seven-factor fitted model.

(b) An underspecified five-factor fitted model.

Figure 8: Rejection rates as well as test set classification accuracies for
approximate and exact C2STs.

C2ST results for the underspecified fitted model in setting (2a) are shown in
Figure 8(b). NN-based C2STs performed well, with classification accuracies
exceeding $0.5$ and rejection rates tending to one as
$N_{\operatorname{test}}$ increases. KNN-based C2STs, however, performed
relatively poorly, with classification accuracies only marginally exceeding
$0.5$, exact tests only obtaining moderate rejection rates when
$N_{\operatorname{test}}=$10\,000$$, and approximate tests always failing to
reject $H_{0}$. These results mirror NNs’ superior performance in Section 4.2.
We also demonstrate the viability of PIs by assessing the underspecified
fitted model’s item-level fit when $N=$10\,000$$. Results in Figure 9 show
that NNs clearly flagged items $17$, $18$, $41$, and $48$ as poor fitting with
all other items having PIs near zero. KNNs also flagged these items, although
PIs were much smaller than for NNs.

Figure 9: Permutation importances for the underspecified five-factor fitted
model when $N=$10\,000$$.

Boxplots of C2ST-RFIs for the overspecified SFM in setting (1b) as well as for
the underspecified FFM in setting (2a) are shown in Figures 10(a) and 10(b),
respectively. C2ST-RFIs for the correctly specified models in settings (1a)
and (2b) were nearly identical to C2ST-RFIs for (1b) and are not shown. The
means of the sampling distributions of both KNN- and NN-based C2ST-RFIs appear
to depend on $N$. C2ST-RFIs for both (1b) and (2a) start out larger than one
when $N_{\operatorname{test}}=750$, with RFIs for (1b) tending to one and RFIs
for (2a) tending to values smaller than one as $N_{\operatorname{test}}$
increases. These results suggest that the SFM fit its data almost perfectly
and that the FFM fit its data relatively poorly, although classifiers
overfitted to noise in the training data for small $N_{\operatorname{train}}$.
The provisory cutoff of $0.9$ obtains similar rejection rates to those
obtained by C2STs: essentially no C2ST-RFIs fall below the threshold in
setting (1a) (i.e., the rejection rate remains near zero), all NN-based C2ST-
RFIs in setting (2b) fall below the threshold as $N_{\operatorname{test}}$
increases (i.e., the rejection rate tends to one), and KNN-based C2ST-RFIs in
setting (2b) start falling below the threshold for large
$N_{\operatorname{test}}$ (i.e., the rejection rate slowly increases).

(a) An overspecified seven-factor fitted model. (b) An underspecified five-
factor fitted model.

Figure 10: C2ST-RFIs for simulated data. The provisory cutoff of $0.9$ is
marked with a solid line.

We assessed each GOF assessment method’s computational efficiency by computing
the total time required to sample all synthetic data and to fit (and possibly
tune hyperparameters for) all associated classifiers. Total run times for
C2STs and for C2ST-RFIs computed in simulation setting (2a) are shown in
Figures 11(a) and 11(b), respectively. Run times for other simulation settings
were very similar and are not shown. Run times increase with increasing
$N_{\operatorname{test}}$, with KNNs always running faster than NNs. For
instance, when $N_{\operatorname{test}}=$10\,000$$, KNN-based C2STs and C2ST-
RFIs have median run times around $4$ and $9$ seconds, respectively, whereas
their NN-based counterparts have median run times around $18$ and $40$
seconds. This finding is likely attributable to the fact that no KNN
hyperparameters are tuned during fitting as well as to the fact that $N$ is
never extremely large (e.g., $N\geq$100\,000$$), in which case NNs would
potentially be much faster than KNNs. Additionally, absolute run times are
never prohibitively large, suggesting that the better-performing NNs are still
a viable choice in practice.

(a) C2ST run times. (b) C2ST-RFI run times.

Figure 11: Run times for KNN- and NN-based C2STs and C2ST-RFIs from the
simulation setting with the underspecified five-factor model.

##  5 Discussion

This work is concerned with the theoretical properties and empirical
performance of machine learning-based parameter estimation and goodness-of-fit
assessment methods for large-scale confirmatory item factor analysis. An
importance-weighted amortized variational estimator implemented via a deep
learning algorithm demonstrated improved parameter recovery as the sample size
increased and obtained comparable estimates to those obtained by the MML
estimator implemented via the MH-RM algorithm. The deep learning algorithm’s
computational efficiency appears to be essentially independent of the sample
size, enabling fitting even with extremely large samples. A simulation-based
test of exact fit called the classifier two-sample test was described and
extended into a test of approximate fit as well a relative fit index.
Approximate C2STs and C2ST-RFIs successfully identified when an IFA model was
correctly or incorrectly specified. C2ST results suggested that more flexible
(but potentially harder to tune) classifiers such as NNs performed better than
less flexible (but potentially easier to tune) classifiers such as KNNs. A
permutation importance technique was demonstrated as a tool for exploring
piece-wise model fit.

The proposed methods have a number of limitations and extensions that may be
addressed in future work.

First, optimization for I-WAVE may be improved to ensure both (1) stabler and
(2) faster convergence. Regarding (1), we note that I-WAVE occasionally
converged to poor local maxima, particularly with small sample sizes. We
handled this issue by refitting the model with new random starting values. We
also manually reduced the learning rate optimization hyperparameter when the
latent dimension was large to ensure stable convergence. It may be feasible to
mitigate instability using more computationally intensive procedures for
choosing starting values such as Cai’s (20102) stochastic EM procedure or Y.
Chen, Li Zhang’s (2019) singular value decomposition procedure. In addition,
the automatic procedures for choosing starting values and for tuning
optimization hyperparameters recommended by Agrawal . (2020) in the context of
general amortized importance-weighted VI may help reduce user burden. At
present, we recommend trying multiple random starts to mitigate convergence to
poor local maxima. We note that this is typically considered good practice for
IFA models fitted via most popular estimators with non-convex objective
functions (e.g., the MML estimator) since convergence to the global maximum is
not guaranteed in general. Regarding (2), we note that convergence could
potentially be sped up using a variant of natural gradient descent (e.g.,
Hoffman ., 2013; Tang Ranganath, 2019). This would be similar in spirit to
recent work by Liu (2020) in which computational efficiency for a JML
estimator was enhanced via a Riemannian optimization procedure that exploits
the manifold structure of the parameter space.

Second, computing standard errors (SEs) and handling missing data with I-WAVE
were not discussed. Approximate SEs may be obtained by evaluating the observed
information matrix at the parameter estimates obtained by maximizing the IW-
ELBO, then inverting this matrix block-wise (Hui ., 2017). As noted by Urban
Bauer (2021), SEs will likely be small for the large-scale applications
considered here. Procedures for handling missing-at-random data (Mattei
Frellsen, 2019) as well as missing-not-at-random data (Ipsen ., 2021) have
been developed for general amortized importance-weighted VI and can be
straightforwardly utilized for I-WAVE.

Third, using approximate C2STs to assess person and piece-wise fit may be more
thoroughly investigated. In terms of person fit, the approach suggested in
Sect. 3.1 remains to be investigated via simulation studies. In terms of
piece-wise fit, the permutation importance measures applied here have
shortcomings including (1) unknown sampling distributions and (2) degraded
performance when predictors are highly correlated (e.g., Hooker Mentch, 2019).
For (1), future work may explore whether accurate $p$-values and confidence
intervals for estimated PIs can be obtained via parametric approximations
(Altmann ., 2010) or bootstrap resampling. For (2), it would likely be
fruitful to explore alternative approaches based on Shapley additive
explanations (SHAP; Lundberg Lee, 2017). SHAP is a game theoretic method for
explaining fitted classifiers’ predictions and has several extensions that
help prevent its performance from degrading in the presence of
multicollinearity (Aas ., 2021; Basu Maji, 2020; Sellereite Jullum, 2020). We
note, however, that SHAP is less computationally efficient than PI,
potentially hampering its application to very large-scale data.

Fourth, C2ST-RFIs may be improved by developing better (1) cutoff criteria,
(2) baseline models, and (3) model complexity penalties. For (1), we studied
the performance of C2ST-RFIs under a fairly narrow range of modeling
conditions. Although our provisory cutoff criterion of $0.9$ performed well in
these conditions, further research is needed to develop more robust cutoff
criteria for evaluating model fit in practice. For (2), we focused on zero-
factor or “independence” baseline models to closely parallel how fit is
assessed in the SEM literature. However, Maydeu-Olivares (20132) notes that
zero-factor baseline models may be overly unrealistic in that they disregard
inter-item correlations, potentially inflating our sense of the proposed
model’s relative value. A more realistic alternative is suggested by Urban
Bauer’s (2021, Section 7.4) note that I-WAVE can be used to fit a fully
nonlinear factor model in which the inverse link function between the observed
and latent variables is approximated using an NN: first fit the proposed model
as well as the fully nonlinear baseline model, then construct a “badness-of-
fit” index that measures the proportional decrement in fit obtained by moving
from the highly flexible baseline to the less flexible model of interest
(e.g., using the ratio
$\widehat{\operatorname{\textit{acc}}}_{\operatorname{base}}/\widehat{\operatorname{\textit{acc}}}_{\operatorname{prop}}$).
For (3), we employed a complexity penalty based on the number of fitted
parameters. Although this parameter counting approach is fast, alternative
penalties that take into account IFA models’ functional forms would provide
more accurate characterizations of complexity that could be used to compare
different models with the same number of parameters (e.g., Bonifay Cai, 2017).

In summary, the methods considered in this work may provide feasible and
promising frameworks for testing hypotheses about the latent structure
underlying large-scale item response data in a computationally efficient
manner. Both I-WAVE and C2STs are highly flexible frameworks that may be
extended in a variety of ways. Some of these extensions are discussed above,
some are discussed in the vast machine learning literature, and some are yet
to be conceived. We view this work as part of a dialogue between machine
learning and psychometrics that is leading to the development of new
extensions and applications with the potential to positively impact both
fields.

## References

  * Aas . (2021) Aas2021Aas, K., Jullum, M. Løland, A.  2021\.  Explaining individual predictions when features are dependent: More accurate approximations to Shapley values Explaining individual predictions when features are dependent: More accurate approximations to Shapley values.  Artificial Intelligence2981–28.  10.1016/j.artint.2021.103502
  * Agrawal . (2020) Agrawal2020Agrawal, A., Sheldon, D. Domke, J.  2020\.  Advances in black-box VI: Normalizing flows, importance weighting, and optimization Advances in black-box VI: Normalizing flows, importance weighting, and optimization.  34th Conference on Advances in Neural Information Processing Systems (NeurIPS 2020).  http://arxiv.org/abs/2006.10343
  * Altmann . (2010) Altmann2010Altmann, A., Toloşi, L., Sander, O. Lengauer, T.  2010\.  Permutation importance: A corrected feature importance measure Permutation importance: A corrected feature importance measure.  Bioinformatics26101340–1347.  10.1093/bioinformatics/btq134
  * Anderson Rubin (1956) Anderson1956Anderson, TW. Rubin, H.  1956\.  Statistical inference in factor analysis Statistical inference in factor analysis.  Proceedings of the Third Berkeley Symposium on Mathematical Statistics and ProbabilityV111.
  * Basu Maji (2020) Basu2020Basu, I. Maji, S.  2020\.  Multicollinearity correction and combined feature effect in Shapley values Multicollinearity correction and combined feature effect in Shapley values.  arXiv.  http://arxiv.org/abs/2011.01661
  * Béguin Glas (2001) Beguin2001Béguin, AA. Glas, CA.  2001\.  MCMC estimation and some model-fit analysis of multidimensional IRT models MCMC estimation and some model-fit analysis of multidimensional IRT models.  Psychometrika664541–562.  10.1007/BF02296195
  * Bentler (1990) Bentler1990Bentler, PM.  1990\.  Comparative fit indices in structural equation models Comparative fit indices in structural equation models.  Psychological Bulletin1072238–246.
  * Bentler (1995) Bentler1995Bentler, PM.  1995\.  EQS structural equations program manual. EQS structural equations program manual.  Encino, CA.
  * Bentler Bonett (1980) Bentler1980aBentler, PM. Bonett, DG.  1980\.  Significance tests and goodness of fit in the analysis of covariance structures Significance tests and goodness of fit in the analysis of covariance structures.  Psychological Bulletin883588–606.  10.1037/0033-2909.88.3.588
  * Biernacki . (2003) Biernacki2003Biernacki, C., Celeux, G. Govaert, G.  2003\.  Choosing starting values for the EM algorithm for getting the highest likehood in multivariate Gaussian mixture models Choosing starting values for the EM algorithm for getting the highest likehood in multivariate Gaussian mixture models.  Computational Statistics and Data Analysis413-4561–575.  10.1016/S0167-9473(02)00163-9
  * Blei . (2017) Blei2017Blei, DM., Kucukelbir, A. McAuliffe, JD.  2017\.  Variational inference: A review for statisticians Variational inference: A review for statisticians.  Journal of the American Statistical Association112518859–877.  10.1080/01621459.2017.1285773
  * Bock Aitkin (1981) Bock1981Bock, RD. Aitkin, M.  1981\.  Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm.  Psychometrika464443–459.  10.1007/BF02293801
  * Bock . (1988) Bock1988Bock, RD., Gibbons, R. Muraki, E.  1988\.  Full-information item factor analysis Full-information item factor analysis.  Applied Psychological Measurement123261–280.  10.1177/014662168801200305
  * Bollen (19891) BollenFit1989Bollen, KA.  19891\.  A new incremental fit index for general structural equation models A new incremental fit index for general structural equation models.  Sociological Methods & Research173303–316.
  * Bollen (19892) Bollen1989Bollen, KA.  19892\.  Structural equations with latent variables Structural equations with latent variables.  John Wiley & Sons, Inc.
  * Bonifay Cai (2017) Bonifay2017Bonifay, W. Cai, L.  2017\.  On the complexity of item response theory models On the complexity of item response theory models.  Multivariate Behavioral Research524465–484.  https://doi.org/10.1080/00273171.2017.1309262 10.1080/00273171.2017.1309262
  * Bottou . (2018) Bottou2018Bottou, L., Curtis, FE. Nocedal, J.  2018\.  Optimization methods for large-scale machine learning Optimization methods for large-scale machine learning.  SIAM Review602223–311.  10.1137/16M1080173
  * Breiman (2001) BreimanRF2001Breiman, L.  2001\.  Random forests Random forests.  Machine Learning4515–32.  10.1007/9781441993267_5
  * Burda . (2016) Burda2016Burda, Y., Grosse, R. Salakhutdinov, R.  2016\.  Importance weighted autoencoders Importance weighted autoencoders.  4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings1–14.
  * Butcher . (1989) Butcher1989Butcher, JN., Dahlstrom, W., Graham, J., Tellegen, A. Kaemmer, B.  1989\.  MMPI-2: Manual for administration and scoring MMPI-2: Manual for administration and scoring.  Minneapolis, MNUniversity of Minnesota Press.
  * Cai (20101) Cai2010aCai, L.  20101\.  High-dimensional exploratory item factor analysis by a Metropolis-Hastings Robbins-Monro algorithm High-dimensional exploratory item factor analysis by a Metropolis-Hastings Robbins-Monro algorithm.  Psychometrika75133–57.  10.1007/s11336-009-9136-x
  * Cai (20102) Cai2010Cai, L.  20102\.  Metropolis-Hastings Robbins-Monro algorithm for confirmatory item factor analysis Metropolis-Hastings Robbins-Monro algorithm for confirmatory item factor analysis.  Journal of Educational and Behavioral Statistics353307–335.  10.3102/1076998609353115
  * Cai (2017) Cai2017Cai, L.  2017\.  flexMIRT®: Flexible multilevel multidimensional item analysis and test scoring. flexMIRT®: Flexible multilevel multidimensional item analysis and test scoring.  Chapel Hill, NCVector Psychometric Group.
  * Cai Hansen (2013) Cai2013Cai, L. Hansen, M.  2013\.  Limited-information goodness-of-fit testing of hierarchical item factor models Limited-information goodness-of-fit testing of hierarchical item factor models.  British Journal of Mathematical and Statistical Psychology662245–276.  10.1111/j.2044-8317.2012.02050.x
  * Chalmers (2012) Chalmers2012Chalmers, RP.  2012\.  Mirt: A multidimensional item response theory package for the R environment Mirt: A multidimensional item response theory package for the R environment.  Journal of Statistical Software48April 2012.  10.18637/jss.v048.i06
  * X. Chen . (2019) Chen2019Chen, X., Liu, S., Sun, R. Hong, M.  2019\.  On the convergence of a class of ADAM-type algorithms for non-convex optimization On the convergence of a class of ADAM-type algorithms for non-convex optimization.  Internation Conference on Learning Representations.
  * Y. Chen, Filho . (2019) ChenBeta32019Chen, Y., Filho, TS., Prudêncio, RBC., Diethe, T. Flach, P.  2019\.  $\beta$3-IRT : A new item response model and its applications $\beta$3-IRT : A new item response model and its applications.  Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS 2019)89.
  * Y. Chen, Li Zhang (2019) ChenCJMLE2019Chen, Y., Li, X. Zhang, S.  2019\.  Joint maximum likelihood estimation for high-dimensional exploratory item factor analysis Joint maximum likelihood estimation for high-dimensional exploratory item factor analysis.  Psychometrika841124–146.  10.1007/s11336-018-9646-5
  * Cho . (2020) Cho2020Cho, AE., Wang, C., Zhang, X. Xu, G.  2020\.  Gaussian variational estimation for multidimensional item response theory Gaussian variational estimation for multidimensional item response theory.  British Journal of Mathematical and Statistical Psychology.  10.1111/bmsp.12219
  * Chwialkowski . (2015) Chwialkowski2015Chwialkowski, K., Ramdas, A., Sejdinovic, D. Gretton, A.  2015\.  Fast two-sample testing with analytic representations of probability measures Fast two-sample testing with analytic representations of probability measures.  Advances in Neural Information Processing Systems1981–1989.
  * Cremer . (2018) Cremer2018Cremer, C., Li, X. Duvenaud, D.  2018\.  Inference suboptimality in variational autoencoders Inference suboptimality in variational autoencoders.  35th International Conference on Machine Learning, ICML 201831749–1760.
  * Cudeck Henly (1991) Cudeck1991Cudeck, R. Henly, SJ.  1991\.  Model selection in covariance structures analysis and the ”problem” of sample size: A clarification Model selection in covariance structures analysis and the ”problem” of sample size: A clarification.  Psychological Bulletin1093512–519.  10.1037/0033-2909.109.3.512
  * Curi . (2019) Curi2019Curi, M., Converse, GA. Hajewski, J.  2019\.  Interpretable variational autoencoders for cognitive models Interpretable variational autoencoders for cognitive models.  International Joint Conference on Neural Networks International joint conference on neural networks ( 1–8).  Budapest, Hungary.  10.1109/IJCNN.2019.8852333
  * Cybenko (1989) Cybenko1989Cybenko, G.  1989\.  Approximation by superpositions of a sigmoidal function Approximation by superpositions of a sigmoidal function.  Mathematics of Control, Signals, and Systems2303–314.
  * Domke Sheldon (2018) Domke2018Domke, J. Sheldon, D.  2018\.  Importance weighting and variational inference Importance weighting and variational inference.  Advances in Neural Information Processing Systems2018-Decem4470–4479.
  * Duda . (2001) Duda2001Duda, RO., Hart, PE. Stork, DG.  2001\.  Pattern classification Pattern classification.  New York, NYWiley.
  * Edwards (2010) Edwards2010Edwards, MC.  2010\.  A Markov chain Monte Carlo approach to confirmatory item factor analysis A Markov chain Monte Carlo approach to confirmatory item factor analysis.  Psychometrika753474–497.  10.1007/s11336-010-9161-9
  * Ehm (1991) Ehm1991Ehm, W.  1991\.  Binomial approximation to the Poisson binomial distribution Binomial approximation to the Poisson binomial distribution.  Statistics & Probability Letters117–16.
  * Fisher . (2019) Fisher2019Fisher, A., Rudin, C. Dominici, F.  2019\.  All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously.  Journal of Machine Learning Research20.
  * Fix Hodges (1951) Fix1951Fix, E. Hodges, J.  1951\.  Discriminatory analysis, nonparametric discrimination: Consistency properties Discriminatory analysis, nonparametric discrimination: Consistency properties .  Randolph Field, TXUSAF School of Aviation Medicine.
  * Friedman (2003) Friedman2003Friedman, JH.  2003\.  On multivariate goodness of fit and two sample testing On multivariate goodness of fit and two sample testing.  Proceedings of the Conference on Statistical Problems in Particle Physics, Astrophysics, and Cosmology (PHYSTAT2003)311–313.
  * Ghosh . (2020) Ghosh2020Ghosh, RP., Mallick, B. Pourahmadi, M.  2020\.  Bayesian estimation of correlation matrices of longitudinal data Bayesian estimation of correlation matrices of longitudinal data.  Bayesian Analysis1–20.  10.1214/20-ba1237
  * Glorot Bengio (2010) Glorot2010Glorot, X. Bengio, Y.  2010\.  Understanding the difficulty of training deep feedforward neural networks Understanding the difficulty of training deep feedforward neural networks.  Journal of Machine Learning Research9249–256.
  * Goldberg (1999) Goldberg1999Goldberg, LR.  1999\.  A broad-bandwidth, public domain, personality inventory measuring the lower-level facets of several five-factor models A broad-bandwidth, public domain, personality inventory measuring the lower-level facets of several five-factor models.  I. Meervielde, IJ. Deary, F. De Fruyt F. Ostendorf (), Personality Psychology in Europe Personality psychology in europe ( 7, 7–28).  Tilburg University Press.
  * Goodfellow . (2016) Goodfellow2016Goodfellow, I., Bengio, Y. Courville, A.  2016\.  Deep learning Deep learning.  Cambridge, MAMIT Press.  http://www.deeplearningbook.org
  * Gretton . (2012) Gretton2012Gretton, A., Borgwardt, KM., Rasch, MJ., Schölkopf, B. Smola, A.  2012\.  A kernel two-sample test A kernel two-sample test.  Journal of Machine Learning Research13723–773.
  * Hastie . (2009) Hastie2009Hastie, T., Tibshirani, R. Friedman, J.  2009\.  The elements of statistical learning: Data mining, inference, and prediction The elements of statistical learning: Data mining, inference, and prediction (2nd ed. ).  New YorkSpringer.  https://web.stanford.edu/{~}hastie/Papers/ESLII.pdf
  * Hoffman . (2013) Hoffman2013Hoffman, MD., Blei, DM., Wang, C. Paisley, J.  2013\.  Stochastic variational inference Stochastic variational inference.  Journal of Machine Learning Research141303–1347.
  * Hooker Mentch (2019) Hooker2019Hooker, G. Mentch, L.  2019\.  Please stop permuting features: An explanation and alternatives Please stop permuting features: An explanation and alternatives.  arXiv.  http://arxiv.org/abs/1905.03151
  * Hu Bentler (1999) Hu1999Hu, LT. Bentler, PM.  1999\.  Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives.  Structural Equation Modeling611–55.  10.1080/10705519909540118
  * Huber . (2004) Huber2004Huber, P., Ronchetti, E. Victoria-Feser, MP.  2004\.  Estimation of generalized linear latent variable models Estimation of generalized linear latent variable models.  Journal of the Royal Statistical Society. Series B: Statistical Methodology664893–908.  10.1111/j.1467-9868.2004.05627.x
  * Hui . (2017) Hui2017Hui, FK., Warton, DI., Ormerod, JT., Haapaniemi, V. Taskinen, S.  2017\.  Variational approximations for generalized linear latent variable models Variational approximations for generalized linear latent variable models.  Journal of Computational and Graphical Statistics26135–43.  http://dx.doi.org/10.1080/10618600.2016.1164708 10.1080/10618600.2016.1164708
  * Ipsen . (2021) Ipsen2020Ipsen, NB., Mattei, PA. Frellsen, J.  2021\.  not-MIWAE: Deep generative modelling with missing not at random data not-MIWAE: Deep generative modelling with missing not at random data.  Internation Conference on Learning Representations (ICLR)1–18.  http://arxiv.org/abs/2006.12871
  * Jitkrittum . (2016) Jitkrittum2016Jitkrittum, W., Szabó, Z., Chwialkowski, K. Gretton, A.  2016\.  Interpretable distribution features with maximum testing power Interpretable distribution features with maximum testing power.  Advances in Neural Information Processing Systems181–189.
  * Jöreskog Moustaki (2001) Joreskog2001Jöreskog, KG. Moustaki, I.  2001\.  Factor analysis of ordinal variables: A comparison of three approaches Factor analysis of ordinal variables: A comparison of three approaches.  Multivariate Behavioral Research363347–387.  10.1207/S15327906347-387
  * Kim . (2021) Kim2021Kim, I., Ramdas, A., Singh, A. Wasserman, L.  2021\.  Classification accuracy as a proxy for two-sample testing Classification accuracy as a proxy for two-sample testing.  Annals of Statistics491411–434.  10.1214/20-aos1962
  * Kingma Welling (2014) Kingma2014Kingma, DP. Welling, M.  2014\.  Auto-encoding variational Bayes Auto-encoding variational Bayes.  International Conference on Learning Representations1–14.  http://arxiv.org/abs/1312.6114
  * Kolmogorov (1933) Kolmogorov1933Kolmogorov, AN.  1933\.  Sulla determinazione empirica di una lgge di distribuzione. Sulla determinazione empirica di una lgge di distribuzione.  Inst. Ital. Attuari,Giorn.483–91.
  * Kuhfeld (2019) Kuhfeld2019Kuhfeld, M.  2019\.  A posterior predictive model checking method assuming posterior normality for item response theory A posterior predictive model checking method assuming posterior normality for item response theory.  Applied Psychological Measurement432125–142.  10.1177/0146621618779985
  * LeCun . (2015) LeCun2015LeCun, Y., Bengio, Y. Hinton, G.  2015\.  Deep learning Deep learning.  Nature Methods5211436–444.  10.1038/nmeth.3707
  * Lee . (2016) LeePPMC2016Lee, T., Cai, L. Kuhfeld, M.  2016\.  A poor person’s posterior predictive checking of structural equation models A poor person’s posterior predictive checking of structural equation models.  Structural Equation Modeling232206–220.  http://dx.doi.org/10.1080/10705511.2015.1014041 10.1080/10705511.2015.1014041
  * Levy . (2009) Levy2009Levy, R., Mislevy, RJ. Sinharay, S.  2009\.  Posterior predictive model checking for multidimensionality in item response theory Posterior predictive model checking for multidimensionality in item response theory.  Applied Psychological Measurement337519–537.  10.1177/0146621608329504
  * Liu (2020) Liu2020Liu, Y.  2020\.  A Riemannian optimization algorithm for joint maximum likelihood estimation of high-dimensional exploratory item factor analysis A Riemannian optimization algorithm for joint maximum likelihood estimation of high-dimensional exploratory item factor analysis.  Psychometrika852439–468.  https://doi.org/10.1007/s11336-020-09711-8 10.1007/s11336-020-09711-8
  * Lopez-Paz Oquab (2017) Lopez-Paz2017Lopez-Paz, D. Oquab, M.  2017\.  Revisiting classifier two-sample tests Revisiting classifier two-sample tests.  5th International Conference on Learning Representations, ICLR 2017 - Conference Track ProceedingsSection 31–15.
  * Lundberg Lee (2017) Lundberg2017Lundberg, SM. Lee, SI.  2017\.  A unified approach to interpreting model predictions A unified approach to interpreting model predictions.  31st Conference on Neural Information Processing Systems (NIPS 2017)4766–4775.
  * Maccallum Tucker (1991) MacCallum1991Maccallum, RC. Tucker, LR.  1991\.  Representing sources of error in the common-factor model: Implications for theory and practice Representing sources of error in the common-factor model: Implications for theory and practice.  Psychological Bulletin1093502–511.  file:///Users/Njaal/Dropbox/Library.papers3/Articles/2004/Unknown/2004-8.pdf{%}5Cnpapers3://publication/uuid/1DF22C4D-FB53-4CF9-965D-29251B1CB559
  * Mann Whitney (1947) Mann1947Mann, HB. Whitney, DR.  1947\.  On a test of whether one of two random variables is stochastically larger than the other On a test of whether one of two random variables is stochastically larger than the other.  The Annals of Mathematical Statistics18150–60.
  * Mattei Frellsen (2019) Mattei2019Mattei, PA. Frellsen, J.  2019\.  MIWAE: Deep generative modelling and imputation of incomplete data sets MIWAE: Deep generative modelling and imputation of incomplete data sets.  36th International Conference on Machine Learning, ICML 20197762–7772.
  * Maydeu-Olivares (20131) Maydeu-Olivares2013Maydeu-Olivares, A.  20131\.  Goodness-of-fit assessment of item response theory models Goodness-of-fit assessment of item response theory models.  Measurement11371–101.  10.1080/15366367.2013.831680
  * Maydeu-Olivares (20132) Maydeu-OlivaresWhyFit2013Maydeu-Olivares, A.  20132\.  Why should we assess the goodness-of-fit of IRT models? Why should we assess the goodness-of-fit of IRT models?  Measurement113127–137.  10.1080/15366367.2013.841511
  * Maydeu-Olivares Joe (2005) Maydeu-Olivares2005Maydeu-Olivares, A. Joe, H.  2005\.  Limited- And full-information estimation and goodness-of-fit testing in 2n contingency tables: A unified framework Limited- And full-information estimation and goodness-of-fit testing in 2n contingency tables: A unified framework.  Journal of the American Statistical Association1004711009–1020.  10.1198/016214504000002069
  * Maydeu-Olivares Joe (2006) Maydeu-Olivares2006Maydeu-Olivares, A. Joe, H.  2006\.  Limited information goodness-of-fit testing in multidimensional contingency tables Limited information goodness-of-fit testing in multidimensional contingency tables.  Psychometrika714713–732.  10.1007/s11336-005-1295-9
  * Maydeu-Olivares Joe (2014) Maydeu-Olivares2014Maydeu-Olivares, A. Joe, H.  2014\.  Assessing approximate fit in categorical data analysis Assessing approximate fit in categorical data analysis.  Multivariate Behavioral Research494305–328.  10.1080/00273171.2014.911075
  * McKinley Reckase (1983) McKinley1983McKinley, R. Reckase, M.  1983\.  An extension of the two-parameter logistic model to the multidimensional latent space An extension of the two-parameter logistic model to the multidimensional latent space .  Iowa City, IAThe American College Testing Program.
  * Meng Schilling (1996) Meng1996Meng, XL. Schilling, S.  1996\.  Fitting full-information item factor models and an empirical investigation of bridge sampling Fitting full-information item factor models and an empirical investigation of bridge sampling.  Journal of the American Statistical Association914351254–1267.  10.1080/01621459.1996.10476995
  * Muthén (1978) Muthen1978Muthén, B.  1978\.  Contributions to factor analysis of dichotomous variables Contributions to factor analysis of dichotomous variables.  Psychometrika434551–560.
  * Muthén (1984) Muthen1984Muthén, B.  1984\.  A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators.  Psychometrika491115–132.
  * Natesan . (2016) Natesan2016Natesan, P., Nandakumar, R., Minka, T. Rubright, JD.  2016\.  Bayesian prior choice in IRT estimation using MCMC and variational Bayes Bayesian prior choice in IRT estimation using MCMC and variational Bayes.  Frontiers in Psychology7SEP1–11.  10.3389/fpsyg.2016.01422
  * Pardos (2017) Pardos2017Pardos, ZA.  2017\.  Big data in education and the models that love them Big data in education and the models that love them.  Current Opinion in Behavioral Sciences18107–113.  10.1016/j.cobeha.2017.11.006
  * Paszke . (2017) Paszke2017Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z.Lerer, A.  2017\.  Automatic differentiation in PyTorch Automatic differentiation in PyTorch.  Workshop on Neural Information Processing Systems.  10.1145/24680.24681
  * Pedregosa . (2011) Pedregosa2011Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O.Duchesnay, E.  2011\.  Scikit-learn: Machine learning in Python Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research122825–2830.
  * Pinheiro Bates (1996) Pinheiro1996Pinheiro, JC. Bates, DM.  1996\.  Unconstrained parametrizations for variance-covariance matrices Unconstrained parametrizations for variance-covariance matrices.  Statistics and Computing63289–296.  10.1007/BF00140873
  * Rabe-Hesketh . (2005) Rabe-Hesketh2005Rabe-Hesketh, S., Skrondal, A. Pickles, A.  2005\.  Maximum likelihood estimation of limited and discrete dependent variable models with nested random effects Maximum likelihood estimation of limited and discrete dependent variable models with nested random effects.  Journal of Econometrics1282301–323.  10.1016/j.jeconom.2004.08.017
  * Rapisarda . (2007) Rapisarda2007Rapisarda, F., Brigo, D. Mercurio, F.  2007\.  Parameterizing correlations: A geometric interpretation Parameterizing correlations: A geometric interpretation.  IMA Journal of Management Mathematics18155–73.  10.1093/imaman/dpl010
  * Reddi . (2018) Reddi2018Reddi, SJ., Kale, S. Kumar, S.  2018\.  On the convergence of ADAM and beyond On the convergence of ADAM and beyond.  International Conference on Learning Representations.  http://arxiv.org/abs/1904.09237
  * Rezende . (2014) Rezende2014Rezende, DJ., Mohamed, S. Wierstra, D.  2014\.  Stochastic backpropagation and approximate inference in deep generative models Stochastic backpropagation and approximate inference in deep generative models.  Proceedings of the 31st International Conference on Machine Learning3221278–1286.
  * Samejima (1969) Samejima1969Samejima, F.  1969\.  Estimation of latent ability using a response pattern of graded scores Estimation of latent ability using a response pattern of graded scores.  Psychometrika351139.  10.1007/BF02290599
  * Schilling Bock (2005) Schilling2005Schilling, S. Bock, RD.  2005\.  High-dimensional maximum marginal likelihood item factor analysis by adaptive quadrature High-dimensional maximum marginal likelihood item factor analysis by adaptive quadrature.  Psychometrika703533–555.  10.1007/s11336-003-1141-x
  * Sellereite Jullum (2020) Sellereite2020Sellereite, N. Jullum, M.  2020\.  shapr: An R-package for explaining machine learning models with dependence-aware Shapley values shapr: An R-package for explaining machine learning models with dependence-aware Shapley values.  Journal of Open Source Software5462027.  10.21105/joss.02027
  * Sinharay . (2006) Sinharay2006Sinharay, S., Johnson, MS. Stern, HS.  2006\.  Posterior predictive assessment of item response theory models Posterior predictive assessment of item response theory models.  Applied Psychological Measurement304298–321.  10.1177/0146621605285517
  * Smirnov (1939) Smirnov1939Smirnov, NV.  1939\.  On the estimation of the discrepancy between empirical curves ofdistribution for two independent samples On the estimation of the discrepancy between empirical curves ofdistribution for two independent samples.  Bull. Math. Univ. Moscou223–14.
  * Song Lee (2005) Song2005Song, XY. Lee, SY.  2005\.  A multivariate probit latent variable model for analyzing dichotomous responses A multivariate probit latent variable model for analyzing dichotomous responses.  Statistica Sinica153645–664.
  * Steiger (1980) Steiger1980Steiger, JH.  1980\.  Statistically based tests for the number of common factors Statistically based tests for the number of common factors.  The Annual Meeting of the Psychometric Society. The annual meeting of the psychometric society.  Iowa City, IA.
  * Student (1908) Student1908Student.  1908\.  The probable error of a mean The probable error of a mean.  Biometrika611–25.
  * Tang Ranganath (2019) Tang2019Tang, D. Ranganath, R.  2019\.  The variational predictive natural gradient The variational predictive natural gradient.  36th International Conference on Machine Learning, ICML 20192019-June10735–10744.
  * Thissen Steinberg (1997) Thissen1997Thissen, D. Steinberg, L.  1997\.  A response model for multiple-choice items A response model for multiple-choice items.  WJ. van der Linden RK. Hambleton (), Handbook of modern item response theory Handbook of modern item response theory ( 51–66).  New York, NYSpringer-Verlag.
  * Tsay Pourahmadi (2017) Tsay2017Tsay, RS. Pourahmadi, M.  2017\.  Modelling structured correlation matrices Modelling structured correlation matrices.  Biometrika1041237–242.  10.1093/biomet/asw061
  * G. Tucker . (2019) Tucker2019Tucker, G., Maddison, CJ., Lawson, D. Gu, S.  2019\.  Doubly reparameterized gradient estimators for Monte Carlo objectives Doubly reparameterized gradient estimators for Monte Carlo objectives.  7th International Conference on Learning Representations, ICLR 201920151–14.
  * LR. Tucker Lewis (1973) Tucker1973Tucker, LR. Lewis, C.  1973\.  A reliability coefficient for maximum likelihood factor analysis A reliability coefficient for maximum likelihood factor analysis.  Psychometrika3811–10.
  * Urban Bauer (2021) UrbanBauer2021Urban, CJ. Bauer, DJ.  2021\.  A deep learning algorithm for high-dimensional exploratory item factor analysis A deep learning algorithm for high-dimensional exploratory item factor analysis.  Psychometrika8611–29.  https://doi.org/10.1007/s11336-021-09748-3 10.1007/s11336-021-09748-3
  * Wilcoxon (1945) Wilcoxon1945Wilcoxon, F.  1945\.  Individual comparisons of grouped data by ranking methods Individual comparisons of grouped data by ranking methods.  Biometrics Bulletin1680–83.  10.1093/jee/39.2.269
  * Wirth Edwards (2007) Wirth2007Wirth, RJ. Edwards, MC.  2007\.  Item factor analysis: Current approaches and future directions Item factor analysis: Current approaches and future directions.  Psychological Methods12158–79.  10.1037/1082-989X.12.1.58
  * Woo . (2020) Woo2020Woo, SE., Tay, L. Proctor, RW. ().  2020\.  Big data in psychological research Big data in psychological research.  American Psychological Association.  https://doi.org/10.1037/0000193-000
  * Wu . (2020) Wu2020Wu, M., Davis, RL., Domingue, BW., Piech, C. Goodman, N.  2020\.  Variational item response theory : Fast, accurate, and expressive Variational item response theory : Fast, accurate, and expressive.
  * Yarkoni Westfall (2017) Yarkoni2017Yarkoni, T. Westfall, J.  2017\.  Choosing prediction over explanation in psychology: Lessons from machine learning Choosing prediction over explanation in psychology: Lessons from machine learning.  Perspectives on Psychological Science1261100–1122.  10.1177/1745691617693393
  * Yun . (2020) Yun2020Yun, J., Lozano, AC. Yang, E.  2020\.  A general family of stochastic proximal gradient methods for deep learning A general family of stochastic proximal gradient methods for deep learning.  arXiv.  http://arxiv.org/abs/2007.07484
  * C. Zhang . (2019) Zhang2019Zhang, C., Butepage, J., Kjellstrom, H. Mandt, S.  2019\.  Advances in variational inference Advances in variational inference.  IEEE Transactions on Pattern Analysis and Machine Intelligence4182008–2026.  10.1109/TPAMI.2018.2889774
  * S. Zhang . (2020) Zhang2020Zhang, S., Chen, Y. Liu, Y.  2020\.  An improved stochastic EM algorithm for large-scale full-information item factor analysis An improved stochastic EM algorithm for large-scale full-information item factor analysis.  British Journal of Mathematical and Statistical Psychology73144–71.  10.1111/bmsp.12153

## Appendix A Proof of Theorem 1

At significance level $\alpha$, the decision threshold for