A Note on Improving Variational Estimation for Multidimensional Item Response Theory

Abstract

Survey instruments and assessments are frequently used in many domains of social science. When the constructs that these assessments try to measure become multifaceted, multidimensional item response theory (MIRT) provides a unified framework and convenient statistical tool for item analysis, calibration, and scoring. However, the computational challenge of estimating MIRT models prohibits its wide use because many of the extant methods can hardly provide results in a realistic time frame when the number of dimensions, sample size, and test length are large. Instead, variational estimation methods, such as Gaussian variational expectation–maximization (GVEM) algorithm, have been recently proposed to solve the estimation challenge by providing a fast and accurate solution. However, results have shown that variational estimation methods may produce some bias on discrimination parameters during confirmatory model estimation, and this note proposes an importance-weighted version of GVEM (i.e., IW-GVEM) to correct for such bias under MIRT models. We also use the adaptive moment estimation method to update the learning rate for gradient descent automatically. Our simulations show that IW-GVEM can effectively correct bias with modest increase of computation time, compared with GVEM. The proposed method may also shed light on improving the variational estimation for other psychometrics models.

Developing, refining, and validating survey questionnaires that measure target latent traits such as personality or cognitive abilities has always been a core agenda in education and psychology, and this focus is also extended to health measurement and culminates in a multi-decade initiative on patient-reported outcome measures. Psychometric methods and tools are an integral part of achieving this focus. When the constructs that these assessments try to measure become increasingly complex, multidimensional item response theory (MIRT), also known as item factor analysis, provides a unified framework and convenient statistical tool for item analysis, calibration, and scoring. However, the increasing scale and complexity of survey designs, especially in large-scale assessments (LSA), require MIRT models with many latent factors. For instance, the English Language Proficiency Assessment for the 21st Century (ELPA21) across two gradebands consists of eight domain-level traits measured by more than 600 items (CRESST, 2017). The existing computational algorithms for fitting high-dimensional MIRT models are insufficient to navigate the massive amount of assessment data, reflected by excessively long computation time and unstable estimation results.

MIRT provides a powerful tool for enriching the information gained in educational assessment (Hartig & Höhler, 2009). For instance, cognitive instructional psychology considers “science knowledge” and “mathematical ability” as highly differentiated theoretical constructs that consist of both basic facts and skills as well as deeper or higher-order understanding (Hamilton et al., 1995; Kupermintz et al., 1995). As another example, the 2003 assessment framework of PISA (OECD, 2003) contains a hierarchy of ability dimensions with general “knowledge and skills” at the highest level, followed by reading, math, science, and problem solving. Then at the lowest level are the sub-domains such as “space and shape,” “change and relationships,” and “quantity” nested within math. Hence, dimensions on different levels vary in their degree of generality and abstraction. Oftentimes, the highest level represents a broad competency level, whereas lower levels represent narrower and more specific abilities. If the intention is to model both the overall and lower-level abilities simultaneously, the model will be high dimensional (Briggs & Wilson, 2003).

Even though the research and development in statistics and psychometrics have provided increasingly sophisticated measurement models to better assess constructs in social sciences, the practice still lags behind (Cai & Hansen, 2018). Unidimensional IRT models continue to dominate the current applications in many domains. One reason is that when the number of items, sample size, and the number of dimensions are all large, the current computational algorithms for MIRT estimation may not be powerful enough to produce results in a reasonable time frame (or ever) (CRESST, 2017). For instance, due to the large number of students and items within each gradeband, the operational analysis approach used for ELPA21 is a two-step approach: in the first step, a unidimensional IRT model is fitted to the item response data for each domain subtest to obtain item parameter estimates; then in the second step, a restricted hierarchical model [i.e., testlet model, Wainer et al. (2007), Gibbons and Hedeker (1992), Cai et al. (2011)] is fitted to estimate the correlations between the four domains (Thissen, 2013). Such a two-step process has two limitations: (1) the item parameter calibration errors are ignored in the second step, and (2) the restricted hierarchical model is only an approximation to the independent-cluster MIRT model. Various full-information methods have been proposed to deal with the computational challenge, which are listed below with pros and cons. The list is by no means exhaustive, but it includes some of the most popular methods that are available in commercial software packages or R packages.1

1. 1. Adaptive Gaussian quadrature. Compared to the regular Gauss-Hermite quadrature [e.g., Bock and Aitkin (1981)], even though the number of quadrature points per dimension is reduced, the total number of quadrature points still increases exponentially with the number of dimensions. Moreover, an extra step is needed to compute the posterior mode and variance of latent factors in each iteration, which adds additional computation costs (Pinheiro & Bates, 1995).

2. 2. Monte Carlo techniques. This family of methods include, for instance, the Monte Carlo EM algorithm (McCulloch, 1997; Wang & Xu, 2015), stochastic EM algorithm (von Davier & Sinharay, 2010; Zhang et al., 2020), or Metropolis–Hastings Robbins–Monro algorithm (Cai 2010a, b). These methods circumvent intractable integrations by sampling from the posterior distributions; however, they may still computationally intensive for complicated high-dimensional models. Fully Bayesian estimation methods, such as Markov chain Monte Carlo (MCMC; Albert, 1992; Patz & Junker, 1999) can also be considered in this category. The Bayesian approach is also computationally costly as it needs a long chain to converge for complex models, though it is preferable with smaller sample sizes.

3. 3. Analytic dimension reduction. For models assuming certain conditional independence among factors (such as the bi-factor models), the conditional independence relations can be used to partition the joint space of all latent variables into smaller subsets. As a result, brute force numerical integration over the joint latent space can be replaced by a sequence of integrations over smaller subsets of latent variables, which helps reduce the computation burden dramatically. This strategy to deal with high-dimensional integration challenges is known as analytic dimension reduction  (Cai et al., 2011; Gibbons & Hedeker, 1992; Rijmen et al., 2008). One limitation, though, is that the algebraic manipulations of the likelihood of a specific model might become very complicated, and they differ for different models (e.g., Cai et al., 2011; Gibbons & Hedeker, 1992). Hence, there is no universal rule that applies to any model.

4. 4. Laplace approximation. This method is based on second-order Taylor expansion of the log-integrand around its mode (Lindstrom & Bates, 1988) such that the high-dimensional integral becomes tractable. This method is a classical and popularly used method for generalized linear mixed-effects models (GLMM), and it is available in many software packages, such as the “lem4” R package (Bates et al., 2014). However, this approximation may not be accurate when the dimension increases to 3 or higher, the sample size is small (Jeon et al., 2017), or the likelihood function is skewed.

Besides the full-information methods above, a recent constraint joint maximum likelihood estimation (CJMLE) was proposed by Chen et al. (2019), which is more computationally efficient than many marginal maximum likelihood methods, and the estimator has the theoretical guarantee to be consistent under high-dimensional settings. Extending CJMLE, the singular value decomposition (SVD)-based estimator was proposed by Zhang et al. (2020a), which further improves the performance of CJMLE. These joint maximum likelihood methods enjoy the low computational cost but sacrifice the flexibility of latent factors by treating them as fixed effects. For instance, it would be hard conceptually to generalize the algorithm to a multiple-group condition in which unbiased estimation of group-specific population distributions is often needed than estimation of individual person’s latent trait as a fixed effect.

In light of the limitations of the above-mentioned methods, variational estimation methods that leverage advances in statistical and machine learning have recently gained increasing interests in psychometrics (Cho et al., 2021, 2022; Jeon et al., 2017). Among numerous variational estimation methods, Rijmen & Jeon (2013) was one of the first to use a variational estimation technique for MIRT models that approximates the likelihood function by a computationally tractable lower bound, but it only studied MIRT models with discrete latent factors. Later, a wide range of studies on variational methods were conducted for the estimation of more complex models  (Hui et al., 2017; Natesan et al., 2016). Recently, Jeon et al. (2017) proposed variational maximization–maximization (VMM) algorithm for the generalized linear mixed models (GLMMs), which outperforms Laplace approximation with a small sample size. However, they rely on some iterative numerical algorithms to attain the solutions in each maximization step, resulting in a slow speed in running the algorithm. To further increase computational efficiency, many researchers brought up variational autoencoder (VAE), a deep learning-based variational method to tackle the estimation problems in MIRT models (Curi et al., 2019; Wu et al., 2020). Extending from VAE, the importance-weighted VAE (IW-VAE) is developed and exhibits competitive performances to other estimation methods  (Liu et al., 2022; Urban & Bauer, 2021) at large sample sizes. However, the two IW-VAE methods lack theoretical support for the consistency of estimators. In addition, although they are powerful in handling large-scale data, their performances in small to medium-sample data may not be as well (see Appendix for more details). Cho et al. (2021, 2022) proposed a Gaussian variational expectation–maximization (GVEM) algorithm, which has shown to be computationally fast and produces comparable and sometimes more accurate parameter estimates than the MH-RM algorithm and than the CJMLE method in high-dimensional exploratory item factor analysis models (i.e., M2PL and M3PL in Cho et al., 2021). Moreover, Cho et al. (2021) proved that the estimated parameters from GVEM algorithm are consistent under the high-dimensional setting. However, we found that directly applying the GVEM algorithm in confirmatory MIRT models would generate relatively large bias on discrimination parameters, especially when the correlations among factors are high and the sample size is not large (please see Sect. 2 for the detailed simulation results). Such a bias issue happens commonly to variational estimation for various statistical models (Bishop, 2006).

To correct the bias in the variational algorithms for MIRT models, we propose an importance-weighted GVEM algorithm (denoted as IW-GVEM hereafter) , which is an extension of GVEM algorithm by performing additional steps after GVEM convergence. The primary idea is to use an importance-weighted variational inference technique to create a tighter variational lower bound to the target, otherwise intractable, marginal likelihood. Because the variational lower bound proposed in Cho et al. (2021, 2022) is replaced by a weighted average based on importance sampling (Domke & Sheldon, 2018), the desirable closed-form solution in the M-step is no longer applicable. Instead, we propose to use Adam (Kingma & Ba, 2014), a popular algorithm for first-order gradient-based optimization. This computationally efficient algorithm updates the objective function stochastically based on adaptive estimates of lower-order moments, and it is especially well suited for large data and complex models. Moreover, different from the IW-VAE methods rooted in deep neural network models where substantial theoretical works on the consistency of the estimators remain to be done, our proposed IW-GVEM is a more transparent method that comes with theoretical guarantees under the high-dimensional setting.

In what follows, this note briefly describes the M2PL model and the original GVEM algorithm and then introduces the IW-GVEM algorithm in Sect. 1, followed by a comprehensive simulation study in Sect. 2. We end the paper with discussions and future directions.

1. Methods

1.1. M2PL

Multidimensional 2PL model is one of the most widely used MIRT models in practice (Reckase, 2009). With M2PL, the item response function of the ith individual to the jth item is modeled by

(1) $\begin{array}{c}P(Y_{\mathrm{ij}}=1\mid {\theta }_i)=\frac{exp(a_j^{\top }{\theta }_i-b_j)}{1+exp(a_j^{\top }{\theta }_i-b_j)},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(Y_{ij}=1 \mid \varvec{\theta }_i)= \frac{\exp ({\varvec{a}}_j^\top \varvec{\theta }_i-b_j)}{1+\exp ({\varvec{a}}_j^\top \varvec{\theta }_i-b_j)}, \end{aligned}$$\end{document}

where $Y_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij}$$\end{document} for $i=1,...,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,...,N$$\end{document} and $j=1,...,J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,...,J$$\end{document} is a binary response, $a_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{a}}_j$$\end{document} denotes a K-dimensional vector of item discrimination parameters for item j, and $b_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} specifies the corresponding difficulty level with item difficulty parameter as $b_j/{\Vert a_j\Vert }_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j/\Vert {\varvec{a}}_j\Vert _2$$\end{document} . Following notations in Cho et al. (2021), we use $Y_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Y}}_i$$\end{document} to denote the response vector of the ith subject, and ${\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_i$$\end{document} to denote the latent trait vector of the ith subject. We write $A=({\alpha }_j,j=1,\cdots ,J)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A = (\varvec{\alpha _j}, j=1,\dots , J)$$\end{document} and $B=(b_j,j=1,\cdots ,J)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B = (b_j, j = 1, \dots , J)$$\end{document} . For model identification, oftentimes the means and variances of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} are fixed as zeros and ones, respectively, and the covariance (which is actually correlation) of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} is freely estimated.

1.2. GVEM

Let $\Delta =\left(A,B,\rho \right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Delta }=({{\varvec{A}},{\varvec{B}},\varvec{\rho }})$$\end{document} denote the set of unknown parameters for M2PL, where $\rho$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\rho }$$\end{document} denotes the correlations of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} . As discussed, the population means of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} are fixed at 0, and the population variances are fixed at 1. The correlations among $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} ’s can be freely estimated. Then, the log-marginal likelihood of responses $Y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Y}}$$\end{document} is

(2) $\begin{array}{c}l(\Delta \mid Y)=\sum _{i=1}^{N}logP(Y_i\mid \Delta )=\sum _{i=1}^{N}log\int \prod _{j=1}^{J}P(Y_{\mathrm{ij}}\mid \Delta ,{\theta }_i)\varphi \left({\theta }_i\right)d{\theta }_i,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l(\varvec{\Delta }\mid {\varvec{Y}})= \sum _{i=1}^N \log P({\varvec{Y}}_{i}\mid \varvec{\Delta })=\sum _{i=1}^N\log \int \prod _{j=1}^J P(Y_{ij}\mid \varvec{\Delta },\varvec{\theta }_i) \phi (\varvec{\theta }_i)d\varvec{\theta }_i, \end{aligned}$$\end{document}

where $\varphi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} denotes a K-dimensional Gaussian distribution of $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} with mean 0 and covariance ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Sigma _{\varvec{\theta }}$$\end{document} . It is the potentially high-dimensional integration in Eq. (2) that makes direct maximization of the log-marginal likelihood computationally prohibitive. The log-likelihood of response $Y$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{Y}$$\end{document} has an equivalent form

$\begin{array}{c}l(\Delta \mid Y)=\sum _{i=1}^N{\int }_{{\theta }_i}logP(Y_i\mid \Delta )\times q_i\left({\theta }_i\right)d{\theta }_i,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l(\varvec{\Delta } \mid {\varvec{Y}})=\sum _{i=1}^N \int _{\varvec{\theta }_i} \log P({\varvec{Y}}_{i}\mid \varvec{\Delta }) \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i, \end{aligned}$$\end{document}

where $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i (\varvec{\theta }_i)$$\end{document} can be any probability density function satisfying ${\int }_{{\theta }_i}{q}_i\left({\theta }_i\right)d{\theta }_i=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _{\varvec{\theta }_i} q_i (\varvec{\theta }_i) d \varvec{\theta }_i = 1$$\end{document} .

The main idea behind variational inference is to approximate the intractable integral in Eq. (2) with a computationally feasible form, known as the evidence lower bound (ELBO; Blei et al., 2017; Ormerod & Wand, 2010). Because $P(Y_i\mid \Delta )=P(Y_i,{\theta }_i\mid \Delta )/P({\theta }_i\mid Y_i,\Delta )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{Y}}_{i}\mid \varvec{\Delta }) = P({\varvec{Y}}_{i}, \varvec{\theta }_i \mid \varvec{\Delta })/ P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })$$\end{document} , we write $l(\Delta \mid Y)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l(\varvec{\Delta } \mid {\varvec{Y}})$$\end{document} as

$\begin{array}{cc}l(\Delta \mid Y)=& \sum _{i=1}^N{\int }_{{\theta }_i}log\frac{P(Y_i,{\theta }_i\mid \Delta )}{P({\theta }_i\mid Y_i,\Delta )}\times q_i\left({\theta }_i\right)d{\theta }_i\\ =& \sum _{i=1}^N{\int }_{{\theta }_i}log\frac{P(Y_i,{\theta }_i\mid \Delta )q_i\left({\theta }_i\right)}{P({\theta }_i\mid Y_i,\Delta )q_i\left({\theta }_i\right)}\times q_i\left({\theta }_i\right)d{\theta }_i\\ =& \sum _{i=1}^N{\int }_{{\theta }_i}log\frac{P(Y_i,{\theta }_i\mid \Delta )}{q_i\left({\theta }_i\right)}\times q_i\left({\theta }_i\right)d{\theta }_i+KL\{q_i\left({\theta }_i\right)\mid P({\theta }_i\mid Y_i,\Delta )\},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l(\varvec{\Delta } \mid {\varvec{Y}})= & {} \sum _{i=1}^N \int _{\varvec{\theta }_i} \log \frac{P({\varvec{Y}}_{i}, \varvec{\theta }_i \mid \varvec{\Delta })}{P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })} \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i \\= & {} \sum _{i=1}^N \int _{\varvec{\theta }_i} \log \frac{P({\varvec{Y}}_{i}, \varvec{\theta }_i \mid \varvec{\Delta })q_i (\varvec{\theta }_i)}{P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })q_i (\varvec{\theta }_i)} \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i \\= & {} \sum _{i=1}^N \int _{\varvec{\theta }_i} \log \frac{P({\varvec{Y}}_{i}, \varvec{\theta }_i \mid \varvec{\Delta })}{q_i (\varvec{\theta }_i)} \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i + KL \{q_i (\varvec{\theta }_i) \mid P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta }) \}, \end{aligned}$$\end{document}

where $KL\{q_i\left({\theta }_i\right)\mid P({\theta }_i\mid Y_i,\Delta )\}={\int }_{{\theta }_i}log\frac{q_i\left({\theta }_i\right)}{P({\theta }_i\mid Y_i,\Delta )}\times q_i\left({\theta }_i\right)d{\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$KL \{q_i (\varvec{\theta }_i) \mid P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta }) \} = \int _{\varvec{\theta }_i} \log \frac{q_i (\varvec{\theta }_i) }{P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })} \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i$$\end{document} is nonnegative. This is because

$\begin{array}{cc}-KL\{q_i\left({\theta }_i\right)\mid P({\theta }_i\mid Y_i,\Delta )\}=& {\int }_{{\theta }_i}log\frac{P({\theta }_i\mid Y_i,\Delta )}{q_i\left({\theta }_i\right)}\times q_i\left({\theta }_i\right)d{\theta }_i\\ ⩽& {\int }_{{\theta }_i}\left(\frac{P({\theta }_i\mid Y_i,\Delta )}{q_i\left({\theta }_i\right)},-,1\right)\times q_i\left({\theta }_i\right)d{\theta }_i\\ ⩽& {\int }_{{\theta }_i}P({\theta }_i\mid Y_i,\Delta )d{\theta }_i-{\int }_{{\theta }_i}{q}_i\left({\theta }_i\right)d{\theta }_i\\ =& 1-1=0\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} - KL \{q_i (\varvec{\theta }_i) \mid P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta }) \}= & {} \int _{\varvec{\theta }_i} \log \frac{P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })}{q_i (\varvec{\theta }_i) } \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i \\\leqslant & {} \int _{\varvec{\theta }_i} \left(\frac{P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta })}{q_i (\varvec{\theta }_i) }-1\right) \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i \\\leqslant & {} \int _{\varvec{\theta }_i} P(\varvec{\theta }_i \mid {\varvec{Y}}_{i}, \varvec{\Delta }) d \varvec{\theta }_i - \int _{\varvec{\theta }_i} q_i (\varvec{\theta }_i) d \varvec{\theta }_i \\ {}= & {} 1- 1 = 0 \end{aligned}$$\end{document}

Therefore, we have a lower bound of log-likelihood that

(3) $\begin{array}{cc}l(\Delta \mid Y)⩾& \sum _{i=1}^N{\int }_{{\theta }_i}log\frac{P(Y_i,{\theta }_i\mid \Delta )}{q_i\left({\theta }_i\right)}\times q_i\left({\theta }_i\right)d{\theta }_i\\ =& \sum _{i=1}^N{E}_{q_i\left({\theta }_i\right)}[log\frac{P(Y_i,{\theta }_i\mid \Delta )}{q_i\left({\theta }_i\right)}]=:ELBO,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} l(\varvec{\Delta } \mid {\varvec{Y}})\geqslant & {} \sum _{i=1}^N \int _{\varvec{\theta }_i} \log \frac{P({\varvec{Y}}_{i}, \varvec{\theta }_i \mid \varvec{\Delta })}{q_i (\varvec{\theta }_i)} \times q_i (\varvec{\theta }_i) d \varvec{\theta }_i \nonumber \\= & {} \sum _{i=1}^N E_{q_i(\varvec{\theta }_i)}\bigg [\log \frac{P({\varvec{Y}}_{i},\varvec{\theta }_i \mid \varvec{\Delta } )}{q_i(\varvec{\theta }_i)}\bigg ] =: ELBO, \end{aligned}$$\end{document}

where the last term ${\sum }_{i=1}^N{E}_{q_i\left({\theta }_i\right)}[log\frac{P(Y_i,{\theta }_i\mid \Delta )}{q_i\left({\theta }_i\right)}]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{i=1}^N E_{q_i(\varvec{\theta }_i)}\bigg [\log \frac{P({\varvec{Y}}_{i},\varvec{\theta }_i \mid \varvec{\Delta } )}{q_i(\varvec{\theta }_i)}\bigg ]$$\end{document} is the ELBO for $l\left(\Delta \right|Y)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l(\varvec{\Delta }|{\varvec{Y}})$$\end{document} in Equation (2). Maximizing the log-marginal likelihood is then approximated by maximizing ELBO, and $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} , the variational distribution, needs to be carefully chosen to minimize the gap between the log-marginal likelihood and its ELBO.

The key is to find $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} so that ELBO approximates the marginal likelihood $l\left(\Delta \right|Y)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l(\varvec{\Delta }|{\varvec{Y}})$$\end{document} as close as possible. Note that when $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} is the posterior density of ${\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_i$$\end{document} , i.e., $q_i\left({\theta }_i\right)=P({\theta }_i\mid Y_i,\Delta )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)=P(\varvec{\theta }_i\mid {\varvec{Y}}_{i}, \varvec{\Delta } )$$\end{document} , maximizing ELBO is equivalent to Bock & Aitkin (1981)’s marginal maximum likelihood/expectation–maximization (MML/EM) algorithm. Instead, as the choice of $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} determines the computational cost and success of the algorithm, Cho et al. (2021, 2022) proposed a choice of $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} that satisfied two criteria: (1) it is easy to maximize, and (2) it approximates the true log-marginal likelihood well. Due to the independence of the students’ responses in general IRT models, $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} is selected for each individual separately. Specifically, under M2PL, the joint distribution of ${\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_i$$\end{document} and $Y_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Y}}_i$$\end{document} is,

(4) $\begin{array}{cc}& logP(Y_i,{\theta }_i\mid \alpha ,b,\rho )\\ & =\sum _{j=1}^J\{Y_{\mathrm{ij}}({\alpha }_j^{\top }{\theta }_i-b_j)+log\frac{1}{1+exp({\alpha }_j^{\top }{\theta }_i-b_j)}\}+log{\varphi }_{\theta }\left({\theta }_i\right)\\ & \ge \sum _{j=1}^{J}log\frac{e^{{\xi }_{\mathrm{ij}}}}{1+e^{{\xi }_{\mathrm{ij}}}}+\sum _{j=1}^J{Y}_{\mathrm{ij}}({\alpha }_j^{\top }{\theta }_i-b_j)+\sum _{j=1}^J\frac{b_j-{\alpha }_j^{\top }{\theta }_i-{\xi }_{\mathrm{ij}}}{2}\\ & -\sum _{j=1}^J\eta \left({\xi }_{\mathrm{ij}}\right)\{{(b_j-{\alpha }_j^{\top }{\theta }_i)}^2-{\xi }_{\mathrm{ij}}^2\}+log{\varphi }_{\theta }\left({\theta }_i\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\log P(Y_{i}, \varvec{\theta }_i \mid \varvec{\alpha }, \varvec{b}, { \varvec{\rho }}) \nonumber \\&\quad = \sum _{j=1}^J \Big \{ Y_{ij} (\varvec{\alpha }_j^\top \varvec{\theta }_i - b_j) + \log \frac{1}{1+\exp (\varvec{\alpha }_j^\top \varvec{\theta }_i - b_j)} \Big \} + \log \phi _{\varvec{\theta }}(\varvec{\theta }_i) \nonumber \\&\quad \ge \sum _{j=1}^J \log \frac{e^{\xi _{ij}}}{1+e^{\xi _{ij}}} + \sum _{j=1}^J Y_{ij}(\varvec{\alpha }_j^\top \varvec{\theta }_i - b_j) + \sum _{j=1}^J \frac{b_j -\varvec{\alpha }_j^\top \varvec{\theta }_i - \xi _{ij}}{2} \nonumber \\&\qquad - \sum _{j=1}^J \eta (\xi _{ij})\{(b_j -\varvec{\alpha }_j^\top \varvec{\theta }_i)^2 - \xi _{ij}^2\} + \log \phi _{\varvec{\theta }}(\varvec{\theta }_i) \end{aligned}$$\end{document}

(5) $\begin{array}{cc}& :=l(Y_i,{\theta }_i\mid \alpha ,b,\rho ,{\xi }_{\mathrm{ij}}),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\quad := l(Y_i, \varvec{\theta }_i \mid \varvec{\alpha }, \varvec{b}, { \varvec{\rho }}, \xi _{ij}), \end{aligned}$$\end{document}

where ${\xi }_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ij}$$\end{document} is the variational parameter for the ith subject, which will be updated iteratively in the M-step of GVEM, and $\eta \left({\xi }_{\mathrm{ij}}\right)={\left(2{\xi }_{i,j}\right)}^{-1}[e^{{\xi }_{i,j}}/(1+e^{{\xi }_{i,j}})-1/2]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta (\xi _{ij}) = (2 \xi _{i, j})^{-1}[e^{\xi _{i, j}} /(1+e^{\xi _{i, j}})-1 / 2]$$\end{document} . The derivation is as follows. Because the difficulty of handling the marginal distribution of $P\left(Y_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{Y}}_i)$$\end{document} mostly comes from the logistic sigmoid function, which makes the integration over $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} not a closed form in the E-step. As a result, Cho et al. (2021) used a local variational approximation method (Jordan, 2004). Denote $x_{\mathrm{ij}}=b_j-{\alpha }_j^{\top }{\theta }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ij} = b_j - \varvec{\alpha }_j^\top \varvec{\theta }_i$$\end{document} , the local variational method gives the following variational lower bound for the sigmoid function:

$\begin{array}{cc}\frac{1}{1+exp({\alpha }_j^{\top }{\theta }_i-b_j)}=& \frac{exp\left(x_{\mathrm{ij}}\right)}{1+exp\left(x_{\mathrm{ij}}\right)}\\ =& \underset{{\xi }_{\mathrm{ij}}}{max}\frac{exp\left({\xi }_{\mathrm{ij}}\right)}{1+exp\left({\xi }_{\mathrm{ij}}\right)}exp\left(\frac{x_{\mathrm{ij}}-{\xi }_{\mathrm{ij}}}{2},-,\eta ,\left({\xi }_{\mathrm{ij}}\right),(x_{\mathrm{ij}}^2-{\xi }_{\mathrm{ij}}^2)\right)\\ ⩾& \frac{exp\left({\xi }_{\mathrm{ij}}\right)}{1+exp\left({\xi }_{\mathrm{ij}}\right)}exp\left(\frac{x_{\mathrm{ij}}-{\xi }_{\mathrm{ij}}}{2},-,\eta ,\left({\xi }_{\mathrm{ij}}\right),(x_{\mathrm{ij}}^2-{\xi }_{\mathrm{ij}}^2)\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{1}{1+\exp (\varvec{\alpha }_j^\top \varvec{\theta }_i - b_j)}= & {} \frac{\exp (x_{ij})}{1+\exp (x_{ij}) } \\= & {} \max _{\xi _{ij}} \frac{\exp (\xi _{ij})}{ 1+\exp (\xi _{ij}) } \exp \left\{ \frac{x_{ij} - \xi _{ij} }{2} - \eta (\xi _{ij} ) (x_{ij}^2 - \xi _{ij}^2) \right\} \\\geqslant & {} \frac{\exp (\xi _{ij})}{ 1+\exp (\xi _{ij}) } \exp \left\{ \frac{x_{ij} - \xi _{ij} }{2} - \eta (\xi _{ij} ) (x_{ij}^2 - \xi _{ij}^2) \right\} , \end{aligned}$$\end{document}

and by applying the above lower bound to Eq. (4), we get Eq. (5), which provides a variational lower bound for $logP(Y_i,{\theta }_i\mid \alpha ,b,\rho )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \log P(Y_{i}, \varvec{\theta }_i \mid \varvec{\alpha }, \varvec{b}, { \varvec{\rho }})$$\end{document} .

By variational inference theory, we can show that the variational distributions $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} (for $i=1,...,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,..., N$$\end{document} ) that minimize the distances between the lower bound and the joint distribution follow a Gaussian distribution with closed-form mean and variance, i.e., $q_i\left({\theta }_i\right)\sim N({\theta }_i\mid {\mu }_i,{\Sigma }_i)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q_i(\varvec{\theta }_i) \sim N(\varvec{\theta }_i \mid \mu _i,\Sigma _i )$$\end{document} where the mean parameter of the normal distribution is

(6) $\begin{array}{c}{\mu }_i={\Sigma }_i\times \sum _{j=1}^J\{2\eta \left({\xi }_{i,j}\right)b_j+Y_{\mathrm{ij}}-\frac{1}{2}\}{\alpha }_j\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu _i ={\Sigma _i} \times \sum _{j=1}^J \Big \{ 2 \eta (\xi _{i,j} ) b_j +Y_{ij} - \frac{1}{2} \Big \} \varvec{\alpha }_j \end{aligned}$$\end{document}

and the covariance matrix is

(7) $\begin{array}{c}({\Sigma }_i{)}^{-1}=({\Sigma }_{\theta }{)}^{-1}+2\sum _{j=1}^J\eta \left({\xi }_{i,j}\right){\alpha }_j{\alpha }_j^{\top }.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \big (\Sigma _i \big )^{-1} = \big ({\Sigma _{\varvec{\theta }}} \big )^{-1} + 2\sum _{j=1}^J\eta (\xi _{i,j} ) \varvec{\alpha }_j \varvec{\alpha }_j^\top . \end{aligned}$$\end{document}

In the confirmatory model estimation, we update population covariance matrix ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\varvec{\theta }}$$\end{document} by

(8) $\begin{array}{c}{\Sigma }_{\theta }=\frac{1}{N}\sum _{i=1}^N({\Sigma }_i+{\mu }_i{\mu }_i^{\top }).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\Sigma }_{\varvec{\theta }}=\frac{1}{N}\sum _{i=1}^N (\varvec{\Sigma }_i+\varvec{\mu }_i\varvec{\mu }_i^\top ). \end{aligned}$$\end{document}

But because we need to fix the diagonal elements of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\varvec{\theta }}$$\end{document} during estimation to fix the scale, we propose to rescale ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\varvec{\theta }}$$\end{document} after the M-step converges, i.e.,

$\begin{array}{c}{\Sigma }_{\theta }^{\ast }={\left({\left(\sqrt{\text{diag}\left({\Sigma }_{\theta }\right)}\right)}^{-1}\right)}^{\top }{\Sigma }_{\theta }{\left(\sqrt{\text{diag}\left({\Sigma }_{\theta }\right)}\right)}^{-1},\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\Sigma }_{\varvec{\theta }}^{*}=\left(\left(\sqrt{\textrm{diag}(\varvec{\Sigma }_{\varvec{\theta }})}\right) ^{-1}\right) ^\top \varvec{\Sigma }_{\varvec{\theta }}\left(\sqrt{\textrm{diag}(\varvec{\Sigma }_{\varvec{\theta }})}\right) ^{-1}, \end{aligned}$$\end{document}

and the discrimination parameter needs to be rescaled accordingly, i.e., ${\alpha }_j^{\ast }={\alpha }_j\sqrt{\text{diag}\left({\Sigma }_{\theta }\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }_j^*=\varvec{\alpha }_j\sqrt{\textrm{diag}(\varvec{\Sigma }_{\varvec{\theta }})}$$\end{document} . For the exploratory analysis, ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\varvec{\theta }}$$\end{document} is fixed at an identity matrix during estimation, and a post hoc rotation will then produce proper nonzero correlations. In the following, we assume that the GVEM algorithm has converged and we fix the variational parameter ${\xi }_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ij}$$\end{document} as the final estimates. In other words, we do not update ${\xi }_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ij}$$\end{document} in the later iterative steps and ${\xi }_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ij}$$\end{document} is fixed at the initialization GVEM step in Algorithm 1.

1.3. Importance Sampling

Referring back to the basic idea underlying variational inference, i.e., the ELBO for log-likelihood of response $l(\Delta \mid Y)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ l(\varvec{\Delta } \mid {\varvec{Y}})$$\end{document} in the inequality (3), it can be seen that a tighter lower bound is attained when $R\equiv P(Y_i,{\theta }_i\mid \Delta )/q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\equiv {P({\varvec{Y}}_{i},\varvec{\theta }_i \mid \varvec{\Delta } )}/{q_i(\varvec{\theta }_i)}$$\end{document} around its mean $P(Y_i\mid \Delta )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{Y}}_{i}\mid \varvec{\Delta })$$\end{document} . Therefore, we can consider different random variables with the same mean that are more concentrated. For example, we can draw M i.i.d. samples from $q\left(z\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\varvec{z}})$$\end{document} , and average the estimates as in importance sampling (IS):

(9) $\begin{array}{c}{R}_M=\frac{1}{M}\sum _{m=1}^M{R}_m=\frac{1}{M}\sum _{m=1}^M\frac{p(x,z_m)}{q\left(z_m\right)},z_m\sim q(·).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} R_M = \frac{1}{M}\sum _{m=1}^M R_m = \frac{1}{M} \sum _{m=1}^M \frac{p({\varvec{x}}, {\varvec{z}}_m)}{ q({\varvec{z}}_m)}, \ {\varvec{z}}_m \sim q(\cdot ). \end{aligned}$$\end{document}

This leads to a tighter “importance-weighted ELBO" (IW-ELBO) on $logP\left(x\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log P ({\varvec{x}})$$\end{document} ,

(10) $\begin{array}{c}{\text{IW-ELBO}}_M=E_{q\left(Z\right)}\left(log,\frac{1}{M},\sum _{m=1}^M,\frac{p(z_m,x)}{q\left(z_m\right)}\right):=L_M\left(x\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {IW-ELBO}_M = E_{q({\varvec{Z}})}\left[ \log \frac{1}{M} \sum _{m=1}^M \frac{p({\varvec{z}}_m, {\varvec{x}})}{q({\varvec{z}}_m)} \right] := {\mathcal {L}}_M({\varvec{x}}). \end{aligned}$$\end{document}

It is shown that $L_M\left(x\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M({\varvec{x}})$$\end{document} converge to $logp\left(x\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log p({\varvec{x}})$$\end{document} as M goes to infinity (Burda et al., 2015), which is summarized in the following result.

Proposition 1

For all M, the lower bounds satisfy

$\begin{array}{c}logp\left(x\right)\ge L_{M+1}\ge L_M.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \log p({\varvec{x}}) \ge {\mathcal {L}}_{M+1} \ge {\mathcal {L}}_M. \end{aligned}$$\end{document}

Moreover, if $p(x,z)/q\left(z\right|x)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p({\varvec{x}}, {\varvec{z}}) / q({\varvec{z}} | {\varvec{x}})$$\end{document} is bounded, then $L_M$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M$$\end{document} approaches $logp\left(x\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log p({\varvec{x}})$$\end{document} as M goes to infinity.

Motivated by this result, we use the importance sampling method and calculate the derivatives of $L_M$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M$$\end{document} to further perform gradient-based optimization. Specifically, denote $w_m=p(x,z_m)/q\left(z_m\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_m = p({\varvec{x}}, {\varvec{z}}_m) / q({\varvec{z}}_m) $$\end{document} , then the derivatives of $L_M$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M$$\end{document} with respect to $\theta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }$$\end{document} are

$\begin{array}{cc}{\nabla }_{\theta }{L}_M\left(x\right)& ={\nabla }_{\theta }{E}_{q\left(Z\right)}\left(log,\frac{1}{M},\sum _{m=1}^M,w_m\right)\\ & =E_{q\left(Z\right)}\left({\nabla }_{\theta },log,\frac{1}{M},\sum _{m=1}^M,w_m\right)\\ & =E_{q\left(Z\right)}\left(\sum _{m=1}^M,{\tilde{w}}_m,{\nabla }_{\theta },log,w_m\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nabla _{\varvec{\theta }}{\mathcal {L}}_M ({\varvec{x}})&= \nabla _{\varvec{\theta }} {E}_{q({\varvec{Z}} )}\left[ \log \frac{1}{M} \sum _{m=1}^M w_m \right] \\&= {E}_{q({\varvec{Z}})} \left[ \nabla _{\varvec{\theta }} \log \frac{1}{M} \sum _{m=1}^M w_m \right] \\&= {E}_{q({\varvec{Z}})} \left[ \sum _{m=1}^M {\tilde{w}}_m \nabla _{\varvec{\theta }} \log w_m \right] , \end{aligned}$$\end{document}

where ${\tilde{w}}_m=w_m/{\sum }_{m^{\prime }=1}^M{w}_{m^{\prime }}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{w}}_m = w_m / \sum _{m'=1}^M w_{m'}$$\end{document} and

(11) $\begin{array}{c}{\nabla }_{\theta }log{w}_m={\nabla }_{\theta }logp(x,z_m)-{\nabla }_{\theta }logq\left(z_m\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nabla _{\varvec{\theta }} \log w_m = \nabla _{\varvec{\theta }} \log p({\varvec{x}}, {\varvec{z}}_m) - \nabla _{\varvec{\theta }} \log q({\varvec{z}}_m). \end{aligned}$$\end{document}

1.4. IW-GVEM

The primary idea of IW-GVEM is to replace Eq. (3) with importance-weighted ELBO as in Eq. (10). That is, for each $i=1,\cdots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i = 1,\dots , N$$\end{document} , we draw M samples from $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i)$$\end{document} for S times:

$\begin{array}{c}{\theta }_i^{(s,m)}\sim q_i\left({\theta }_i\right),\text{for}s=1,\cdots ,S,m=1,\cdots ,M.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \varvec{\theta }_i^{(s, m)} \sim q_i(\varvec{\theta }_i), \text { for } s = 1, \dots , S, m = 1, \dots , M. \end{aligned}$$\end{document}

Define $w_i^{(s,m)}=p(Y_i,{\theta }_i^{(s,m)})/q_i\left({\theta }_i^{(s,m)}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_i^{(s, m)} = p(Y_i, \varvec{\theta }_i^{(s,m)}) / q_i(\varvec{\theta }_i^{(s,m)})$$\end{document} , where $p(Y_i,{\theta }_i^{(s,m)})=P(Y_i,{\theta }_i^{(s,m)}\mid \alpha ,b,\rho )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(Y_i, \varvec{\theta }_i^{(s,m)}) = P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }})$$\end{document} as in Eq. (4), and $q_i\left({\theta }_i^{(s,m)}\right)\sim N({\theta }_i^{(s,m)}\mid {\mu }_i,{\Sigma }_i)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i^{(s,m)}) \sim N(\varvec{\theta }_i^{(s,m)} \mid \mu _i, \Sigma _i)$$\end{document} , then $L_M\left(Y\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M ({\varvec{Y}})$$\end{document} can be approximated by

$\begin{array}{c}{L}_M\left(Y\right)\approx \sum _{i=1}^N\left(\frac{1}{S},\sum _{s=1}^S,\left(log,\frac{1}{M},\sum _{m=1}^M,w_i^{(s,m)}\right)\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathcal {L}}_M ({\varvec{Y}}) \approx \sum _{i=1}^N \left(\frac{1}{S}\sum _{s=1}^S\left[ \log \frac{1}{M}\sum _{m=1}^M w_i^{(s, m)} \right] \right) . \end{aligned}$$\end{document}

Note $w_i^{(s,m)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_i^{(s,m)}$$\end{document} is a function of parameters $({\xi }_i,\alpha ,b,\rho )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\xi _i, \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }})$$\end{document} .

To learn parameters, we use a stochastic gradient ascent method, which needs to calculate the gradients of $L_M\left(Y\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}_M({\varvec{Y}})$$\end{document} . Based on Eq. (11), the gradients can be approximated by

$\begin{array}{c}{\nabla }_{\alpha }{L}_M\left(Y\right)\approx \sum _{i=1}^N\left(\frac{1}{S},\sum _{s=1}^S,\sum _{m=1}^M,{\tilde{w}}_i^{(s,m)},{\nabla }_{\alpha },\left(log,P,(Y_i,{\theta }_i^{(s,m)}\mid \alpha ,b,\rho ),-,{\nabla }_{\alpha },log,q_i,({\theta }_i^{(s,m)}\mid Y_i)\right)\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nabla _{\varvec{\alpha }}{\mathcal {L}}_M ({\varvec{Y}}) \approx \sum _{i=1}^N \left(\frac{1}{S}\sum _{s=1}^S \sum _{m=1}^M \tilde{w}_i^{(s,m)} \nabla _{\varvec{\alpha }} \left[ \log P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }}) - \nabla _{\varvec{\alpha }} \log q_i(\varvec{\theta }_i^{(s,m)} \mid Y_i) \right] \right) , \end{aligned}$$\end{document}

where ${\tilde{w}}_i^{(s,m)}=w_i^{(s,m)}/{\sum }_{m^{\prime }=1}^M{w}_i^{(s,m^{\prime })}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{w}_i^{(s,m)} = w_i^{(s,m)} / \sum _{m'=1}^M w_i^{(s,m')}$$\end{document} . Note that $q_i({\theta }_i^{(s,m)}\mid Y_i)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i(\varvec{\theta }_i^{(s,m)} \mid Y_i)$$\end{document} does not depend on the parameters in the current iteration. Therefore, we only need to calculate ${\tilde{w}}_i^{(s,m)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{w}_i^{(s,m)}$$\end{document} and ${\nabla }_{\alpha }P(Y_i,{\theta }_i^{(s,m)}\mid \alpha ,b,\rho )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla _{\varvec{\alpha }} P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }})$$\end{document} . Similarly we can calculate ${\nabla }_b{L}_M\left(Y\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla _{ {\varvec{b}}}{\mathcal {L}}_M ({\varvec{Y}})$$\end{document} and ${\nabla }_{{\Sigma }_{\theta }}{L}_M\left(Y\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla _{\Sigma _{\varvec{\theta }}}{\mathcal {L}}_M ({\varvec{Y}})$$\end{document} . Specifically, we have

(12) $\begin{array}{cc}{\nabla }_{{\alpha }_j}logP\left(Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho \right)=& {\tilde{w}}_i^{(s,m)}{\nabla }_{{\alpha }_j}\left(log,P,(,Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho ,)\right)\\ =& {\tilde{w}}_i^{(s,m)}\left(\left(Y_{\mathrm{ij}},-,1,+,\frac{1}{1+exp({\alpha }_j^{\top }{\theta }_i^{(s,m)}-b_j)}\right),{\theta }_i^{(s,m)}\right),\\ {\nabla }_{b_j}logP\left(Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho \right)=& {\tilde{w}}_i^{(s,m)}{\nabla }_{b_j}\left(log,P,(,Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho ,)\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nabla _{\varvec{\alpha }_j} \log P\left(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}},\varvec{\rho }\right)= & {} \tilde{w}_i^{(s,m)} \nabla _{\varvec{\alpha }_j} \left[ \log P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }}) \right] \nonumber \\= & {} \tilde{w}_i^{(s,m)} \left[ \left(Y_{ij} - 1 + \frac{1}{1 + \exp (\varvec{\alpha }_j^\top \varvec{\theta }_i^{(s,m)} - b_j)}\right) \varvec{\theta }_i^{(s,m)} \right] , \nonumber \\ \nabla _{b_j} \log P\left(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}},\varvec{\rho }\right)= & {} \tilde{w}_i^{(s,m)} \nabla _{b_j} \left[ \log P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }}) \right] \end{aligned}$$\end{document}

(13) $\begin{array}{cc}=& {\tilde{w}}_i^{(s,m)}\left(1,-,Y_{\mathrm{ij}},-,\frac{1}{1+exp({\alpha }_j^{\top }{\theta }_i^{(s,m)}-b_j)}\right),\\ {\nabla }_{{\Sigma }_{\theta }}logP\left(Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho \right)=& {\tilde{w}}_i^{(s,m)}{\nabla }_{{\Sigma }_{\theta }}\left(log,P,(,Y_i,,,{\theta }_i^{(s,m)},\mid ,\alpha ,,,b,,,\rho ,)\right)\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}= & {} \tilde{w}_i^{(s,m)} \left[ 1 - Y_{ij} - \frac{1}{1 + \exp (\varvec{\alpha }_j^\top \varvec{\theta }_i^{(s,m)} - b_j)} \right] ,\nonumber \\ \nabla _{\Sigma _{\varvec{\theta }}} \log P\left(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}},\varvec{\rho }\right)= & {} \tilde{w}_i^{(s,m)} \nabla _{\Sigma _{\varvec{\theta }}} \left[ \log P(Y_i, \varvec{\theta }_i^{(s,m)} \mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }}) \right] \end{aligned}$$\end{document}

(14) $\begin{array}{cc}=& {\tilde{w}}_i^{(s,m)}\left(\frac{1}{2},{\Sigma }_{\theta },-,\frac{1}{2},{\theta }_i^{(s,m)},{\left({\theta }_i^{(s,m)}\right)}^{\top }\right).\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}= & {} \tilde{w}_i^{(s,m)} \left[ \frac{1}{2}\Sigma _{\varvec{\theta }} - \frac{1}{2} \varvec{\theta }_i^{(s,m)} (\varvec{\theta }_i^{(s,m)})^\top \right] . \end{aligned}$$\end{document}

To summarize, in the $(t+1)th$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(t+1)th$$\end{document} iteration, we perform the following:

1. 1. For $i=1,\cdots ,N$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1, \dots , N$$\end{document} , draw M samples from $q_i\left({\theta }_i\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_i (\varvec{\theta }_i)$$\end{document} for S times.

2. 2. Calculate $w_i^{(s,m)}=P(Y_i,{\theta }_i^{(s,m)}\mid \alpha ,b,\rho )/q_i\left({\theta }_i^{(s,m)}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_i^{(s, m)} = P(Y_i, \varvec{\theta }_i^{(s,m)}\mid \varvec{\alpha }, {\varvec{b}}, { \varvec{\rho }}) / q_i(\varvec{\theta }_i^{(s,m)})$$\end{document} and ${\tilde{w}}_i^{(s,m)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{w}_i^{(s,m)} $$\end{document} $=w_i^{(s,m)}/{\sum }_{m^{\prime }=1}^M{w}_i^{(s,m^{\prime })}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= w_i^{(s,m)} / \sum _{m'=1}^M w_i^{(s,m')}$$\end{document} .

3. 3. Calculate the gradients according to Eqs.  (12), (13), and (14).

Algorithm 1 IW-GVEM for M2PL

Proper learning rate scheduling is important in gradient-based algorithms. In this work, we apply the Adaptive moment estimation (Adam) method (Kingma & Ba, 2014), which has been extensively used in deep learning research and applications, to adjust the learning rate in our training process. In Adam, we compute individual adaptive learning rates for each parameter from estimates of the first and second moments of the gradients. Specifically in the tth iteration, we calculate exponential moving averages of the gradient (denoted as $v_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{v}}_t$$\end{document} ) and the squared gradient (denoted as $s_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{s}}_t$$\end{document} ) with exponential decay rates ${\beta }_1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _1$$\end{document} and ${\beta }_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _2$$\end{document} , respectively. The moving averages can be seen as estimates of the first and second moments of the gradients. Then, we correct these biased exponential moving averages by $1-{\beta }_1^t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - \beta _1^t$$\end{document} and $1-{\beta }_2^t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 - \beta _2^t$$\end{document} , respectively, and update parameters using standardized gradients. The concrete steps of generic Adam are provided below, where $g_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{g}}_t$$\end{document} is the gradient (corresponding to that in Eqs. (12), (13) and (14), respectively) in the tth iteration:

1. 1. $v_t={\beta }_1{v}_{t-1}+(1-{\beta }_1)g_t$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{v}}_t = \beta _1 {\varvec{v}}_{t-1} + (1-\beta _1){\varvec{g}}_t$$\end{document} (update biased first moment estimate)

2. 2. $r_t={\beta }_2{r}_{t-1}+(1-{\beta }_2)g_t^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${ {\varvec{r}}_t} = \beta _2 { {\varvec{r}}_{t-1}} + (1-\beta _2){\varvec{g}}_t^2$$\end{document} (update biased second moment estimate)

3. 3. ${\hat{v}}_t=v_t/(1-{\beta }_1^t),{\hat{r}}_t=r_t/(1-{\beta }_2^t)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{v}}}_t = {\varvec{v}}_t / (1 - \beta _1^t),\ { \hat{{\varvec{r}}}_t} = { {\varvec{r}}_t} / (1 - \beta _2^t)$$\end{document} (compute bias-corrected moment estimates)

4. 4. ${\hat{g}}_t=\eta \widehat{v_t}/(\sqrt{{\hat{r}}_t}+ϵ)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{g}}}_t = \eta \hat{{\varvec{v}}_t} /(\sqrt{{ \hat{{\varvec{r}}}_t}} + \epsilon )$$\end{document} , where $\eta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} is learning rate (update the final gradient)

With this, the proposed importance-weighted Gaussian variational EM (IW-GVEM) algorithm is summarized in Algorithm 1. For the choice of hyperparameters, we follow the suggestions in Kingma & Ba (2014) and adopt the default setting that ${\beta }_1=0.9$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _1 = 0.9$$\end{document} and ${\beta }_2=0.999$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _2 = 0.999$$\end{document} . Empirically in our simulation studies, for better convergence performance, we let the learning rate of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} to be $0.1\eta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.1\eta $$\end{document} while the learning rate for $a$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{a}$$\end{document} and $b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{b}$$\end{document} to be $\eta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} , and we search for an optimal learning rate $\eta$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} with the maximum ELBO over a list $\{0.01,0.05,0.1,0.5\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{0.01, 0.05, 0.1, 0.5\}$$\end{document} . Lastly, we set $ϵ=0.001$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon = 0.001$$\end{document} .

In terms of convergence criteria, we evaluate the Euclidean norm of the difference between the estimated parameters of the current step and those of the previous step. When the difference is less than a certain tolerance value, the algorithm is stopped. For our simulation studies, in obtaining the initial model parameter using the GVEM algorithm, we reach convergence at $(l+1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(l+1)$$\end{document} th iteration if $\Vert {\alpha }_{\mathrm{GV}}^{l+1}-{\alpha }_{\mathrm{GV}}^l{\Vert }_2+\Vert b_{\mathrm{GV}}^{t+1}-b_{\mathrm{GV}}^t{\Vert }_2+{\Vert {\Sigma }_{\theta ,GV}^{l+1}-{\Sigma }_{\theta ,GV}^l\Vert }_2⩽0.0001$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \varvec{\alpha }_{GV}^{l+1} - \varvec{\alpha }_{GV}^{l}\Vert _2 + \Vert \varvec{b}_{GV}^{t+1} - \varvec{b}_{GV}^{t}\Vert _2 + \Vert \varvec{\Sigma }_{\theta , GV}^{l+1} - \varvec{\Sigma }_{\theta , GV}^{l}\Vert _2 \leqslant 0.0001$$\end{document} . In IW-GVEM, we reach convergence at $(t+1)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(t+1)$$\end{document} th iteration when $max\{\Vert {\alpha }^{t+1}-{\alpha }^t{\Vert }_2,\Vert b^{t+1}-b^t{\Vert }_2,\Vert {\Sigma }_{\theta }^{t+1}-{\Sigma }_{\theta }^t{\Vert }_2\}⩽0.0001$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\max \{ \Vert \varvec{\alpha }^{t+1} - \varvec{\alpha }^{t}\Vert _2, \Vert \varvec{b}^{t+1} - \varvec{b}^{t}\Vert _2, \Vert \varvec{\Sigma }_{\theta }^{t+1} - \varvec{\Sigma }_{\theta }^{t}\Vert _2 \} \leqslant 0.0001$$\end{document} or the iteration stops when it reaches certain maximum iteration number.

2. Simulation Studies

2.1. Design

We conducted comprehensive simulation studies to evaluate the performance of the proposed method under various manipulated conditions. We follow similar designs as in Cho et al. (2021) and consider different settings: (1) sample size: N = 200 or 500; (2) number of domains: K = 2 or 5; (3) test length: J = 30 if K = 2 or J = 55 if K = 5; (4) both within and between multidimensional structures; (5) factor correlations: low correlation $r\sim \text{unif}(0.1,0.3)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r\sim \text {unif}(0.1, 0.3)$$\end{document} or high correlation $r\sim \text{unif}(0.5,0.7)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r \sim \text {unif}(0.5, 0.7)$$\end{document} ; and (6) confirmatory or exploratory analysis.

Figure. 1 Bias for $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} under confirmatory analysis.

Figure. 2 RMSE for $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} under confirmatory analysis.

Similar to Cho et al. (2021), for the between-item multidimensional structure, we had equal numbers of items loaded on each factor. For the within-item multidimensional structure, when $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} , about one-third of the items were loaded onto the first, or the second, or both factors, respectively. In the cases where $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} , there were about one-third of the items loaded onto one, two, or three factors, respectively. For the model parameters, we simulated the item discrimination parameters ${\alpha }_{j,k}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{j,k}$$\end{document} from uniform distribution on [1, 2], and difficulty parameter $b_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} from the standard normal distribution. We generated the latent traits ${\theta }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\theta }_j$$\end{document} from multivariate normal distribution $N(0,{\Sigma }_{\theta })$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N({\varvec{0}}, \varvec{\Sigma }_{\varvec{\theta }})$$\end{document} , where the diagonal elements of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\varvec{\theta }}$$\end{document} were all 1 and off-diagonal elements were generated from uniform distributions. Specifically, in high-correlation settings, the uniform distribution was set to be $\text{unif}(0.5,0.7)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {unif}(0.5, 0.7)$$\end{document} , whereas in the low-correlation settings we set it to be $\text{unif}(0.1,0.3)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {unif}(0.1, 0.3)$$\end{document} .

For evaluation, we compared the bias and root-mean-squared errors (RMSEs) of model parameters, as well as computation time between GVEM and IW-GVEM. For exploratory analysis, we did a promax rotation after model convergence, and compared the rotated parameters to the true values (Cho et al., 2022). For IW-GVEM, we first ran GVEM algorithm to get initial estimates of model parameters, and then ran several gradient descent steps using importance sampling to correct the bias. To select a proper initial learning rate for the gradient algorithm, we first sampled a set of data aside based on the GVEM estimates. After we got model parameter estimates using importance sampling, we calculated the lower bound as in our objective function based on the previously sampled data set, and chose the learning rate corresponding to the largest lower bound. In the simulation studies, we set S and M to be 10. Our empirical experiments have shown that increasing S and M did not result in significant improvements and 10 was large enough for the simulation settings. The results were averaged over 100 repetitions.

2.2. Results

Figures 1 and 2 present the bias and RMSE of confirmatory M2PL model when $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} . Note that in confirmatory analysis, there are discrimination parameters specified to be zeros. These zero-constrained terms are excluded in the bias and RMSE computation. The two separately colored boxes represent the distribution of respective criteria across 100 replications from IW-GVEM (denoted as “IS” in the figure) and the original GVEM algorithm. As shown, GVEM already performs well by producing close to 0 bias for $b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{b}}$$\end{document} and ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} . It is the discrimination parameter, $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} , that has a non-ignorable bias. The IW-GVEM algorithm effectively corrects such bias on $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} across all conditions without deteriorating the estimation accuracy of other parameters. And because the bias is corrected, the RMSE of $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} is also smaller consistently compared to that from GVEM, whereas again, there is no appreciable difference between IW-GVEM and GVEM in terms of RMSE on $b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{b}}$$\end{document} and ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} . Sorting through the manipulated conditions, it is “within-item” multidimensional structure in combination with high factor correlation tends to yield larger RMSE for both methods and all parameters.

Figure. 3 Bias for $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} under confirmatory analysis.

Figure. 4 RMSE for $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} under confirmatory analysis.

Figure. 5 Bias for $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} under exploratory analysis.

Figure. 6 RMSE for $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} under exploratory analysis.

Figure. 7 Bias for $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} under exploratory analysis.

Figure. 8 RMSE for $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} under exploratory analysis.

Figures 3 and 4 present the bias and RMSE of confirmatory M2PL model when $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} . The trend observed from the $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} condition continues to hold here. That is, IW-GVEM can correct bias on $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} effectively and hence also brings down its RMSE, whereas bias on the other parameters is already close to 0 from both methods and their RMSEs are also comparable. Increasing the number of dimensions certainly makes the model estimation harder to converge, and the estimates are also more variable, especially for $b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{b}}$$\end{document} and ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} , as reflected by wider boxes for those parameters in Fig. 3.

Figures 5, 6, 7, and 8 presents the results from exploratory estimation condition, in the same order as before. For the exploratory M2PL model estimation, GVEM generally performs well and the bias on $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} is already small to begin with. This is consistent with the results reported in literature (Cho et al., 2021, 2022). Even so, under all settings, the RMSEs of IW-GVEM are still smaller than or equal to that of GVEM. IW-GVEM can still further bring down the bias of $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} to near 0 for most cases. The exceptional case when the bias of $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} from IW-GVEM is larger than the bias from GVEM is for the “within item, correlation is high" condition. This case is the most difficult case where the items were loaded on factors via a more complicated setting and the correlations among factors are relatively high. Nonetheless, this special case has overall good estimation performance as the estimation bias from IW-GVEM is still close to the bias from GVEM, and the RMSE from IW-GVEM is lower than the RMSE from GVEM. In addition, when $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} the bias of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} appears to depart from 0 and IW-GVEM does not correct for such bias, although the RMSE of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} is kept small across the board. The bias of ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} gets closer to 0 when K increases and when the factor correlation is low. Because in the exploratory estimation mode, specific types of rotations will affect resulting factor correlations, the bias in ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\Sigma }_{\theta }$$\end{document} estimation is less of a concern. Although the increase in the number of dimensions K could lead to a more complicated model and bring challenges to parameter estimation, the increase in test length, on the other hand, improves the estimation accuracy of parameters. Specifically, at $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 5$$\end{document} , we use test length $J=55$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J = 55$$\end{document} which is greater than $J=30$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J = 30$$\end{document} at $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 2$$\end{document} . This increase in test length explains the results that the biases at $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K =5$$\end{document} are closer to 0 than that at $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K =2$$\end{document} for some cases. Overall, the results from GVEM and IW-GVEM are very close.

Table 1 presents the computation time for confirmatory M2PL estimation under both GVEM and IW-GVEM algorithms. Understandably, IW-GVEM takes longer time under all conditions because both the important sampling step and the gradient descent optimization are time-consuming compared to closed-form updates in GVEM. Unsurprisingly, both methods need longer time for larger sample sizes. It is more interesting to note that, other things being equal, when the multidimensional structure is “within-item,” GVEM almost doubles (when $K=2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=2$$\end{document} ) or sometimes even triples (when $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} ) the computation time compared to the “between-item” condition. But for IW-GVEM, the computation time is rather stable across these two multidimensional structures. Similarly, high correlation among factors is known to be more challenging, hence computation time increases by about 50% or more for GVEM from low to high-correlation conditions, but the computation time of IW-GVEM seems to be unaffected. These all suggest that IW-GVEM is better suited for more complex models. The same patterns remain for the exploratory M2PL estimation, as shown in Table 2, although exploratory analysis in general takes longer time than confirmatory analysis, simply because more parameters are needed to be updated simultaneously.

Table 1 Computation time (seconds) for the confirmatory M2PL estimation.

Table 2 Computation time (seconds) for the exploratory M2PL estimation.

3. Discussion

In this note, we proposed an importance-weighted version of GVEM to correct its bias on the $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} estimates in the confirmatory M2PL models. Because the evidence lower bound (ELBO), a key component of variational inference, is derived based on Jensen’s inequality, the ELBO will approximate the log-marginal distribution (i.e., $logP\left(X\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log P({\varvec{X}})$$\end{document} ) more closely when $R\equiv P(X,Z)/q\left(Z\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\equiv {P({\varvec{X}},{\varvec{Z}})}/{q({\varvec{Z}})}$$\end{document} is more concentrated around its mean $P\left(X\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{X}})$$\end{document} . Hence, the primary idea of IW-GVEM is to replace R with its sample mean by drawing i.i.d. samples from variational distribution $q\left(z\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q({\varvec{z}})$$\end{document} . In so doing, we achieve a tighter bound of Jensen’s inequality, but at the slight cost of computational efficiency. The added computation time is mainly due to sampling in the E-step and gradient descent in the M-step. From our simulation results, the bias correction is effective for confirmatory models and the extra computation time is acceptable because even with additional computational cost, the total time is still short. In fact, the time increase from GVEM to IW-GVEM is at a slow rate in that the time ratio between the two methods is smaller for more complex models (i.e., $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} , within-item multidimensional structure, and high correlations). Note that for exploratory M2PL models, the original GVEM is still recommended because it already produces almost unbiased results and hence importance sampling seems unnecessary, although it does not introduce any undesirable bias either. Theoretically, Cho et al. (2021) proved that the estimated factor loading matrix and estimated latent factor from the GVEM algorithm is consistent as $N\to \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \rightarrow \infty $$\end{document} and $J\to \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J \rightarrow \infty $$\end{document} . The proposed IW-GVEM algorithm is based on the GVEM estimation, hence with consistent initial GVEM estimators, the final estimators from the IW-GVEM algorithm also have the theoretical guarantee to be consistent in the high-dimensional setting. Moreover, compared to ELBO in GVEM, the importance-weighted ELBOs are greatly improved after importance sampling. In finite-sample simulations, importance-weighted ELBOs at $M=5,10,50$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M = 5, 10, 50$$\end{document} , and 100 are all larger than ELBO from GVEM and converge as M increases (See Appendix B).

In IW-GVEM, we propose to use the adaptive moment estimation method to automatically update the learning rate on the fly. Our preliminary results showed that the Adam algorithm performs better than fixed learning rate. Further, we also evaluated the effect of Monte Carlo sample size (i.e., $S=10,50,100$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=10, 50, 100$$\end{document} ) and sample size for the importance sampling step (i.e., $M=10,50$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M=10, 50$$\end{document} ) and noted essentially the same results. Hence, we set $S=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S=10$$\end{document} and $M=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M=10$$\end{document} in our simulation study, which explains the only modest increase in computation time.

Aside from GVEM, another recently proposed fast algorithm for high-dimensional IRT estimation is the joint maximum likelihood estimation (Chen et al., 2019). This method treats the latent abilities as fixed effect parameters instead of random variables. Although this approach is innovative and their algorithm appears to produce accurate parameter estimates efficiently, the interpretation of person parameters is different such that caution needs to be exercised when one intends to generalize findings to a certain population. Plus, treating each individual as a separate fixed effect is, at the conceptual level, hard to justify when generalizing M2PL to a multiple-group MIRT model. This is because the goal of a multiple-group extension is to allow for unbiased marginal estimation of group-specific population distributions.

Instead, the GVEM method can be generalized to multiple-group MIRT in a more straightforward fashion. Our other study exploring multiple-group GVEM for differential item functioning detection (DIF) reveals that it can very well detect uniform DIF, but the power of detecting DIF on discrimination parameter is low. This is likely due to the estimation bias on $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha }$$\end{document} from GVEM in the confirmatory model estimation, and hence the IW-GVEM will likely improve detection of the non-uniform DIF, in particular the DIF on discrimination parameters. Our study can also be extended in other directions. For instance, like in Cho et al. (2021), the IW-GVEM can be extended to M3PL models. Moreover, the current IW-GVEM algorithm does not automatically output standard error of item parameter estimates, and hence future studies may consider combining it with the supplemented EM algorithm (Cai, 2008; Chen & Wang, 2021) to produce accurate SE estimates. In addition to MIRT, the proposed method may also shed light on improving the performance of the variational estimation for other psychometric models, such as generalized linear mixed models (Jeon et al., 2017) and cognitive diagnosis models (Yamaguchi & Okada, 2020, 2020a).