Variational Estimation for Multidimensional Generalized Partial Credit Model

Abstract

Multidimensional item response theory (MIRT) models have generated increasing interest in the psychometrics literature. Efficient approaches for estimating MIRT models with dichotomous responses have been developed, but constructing an equally efficient and robust algorithm for polytomous models has received limited attention. To address this gap, this paper presents a novel Gaussian variational estimation algorithm for the multidimensional generalized partial credit model. The proposed algorithm demonstrates both fast and accurate performance, as illustrated through a series of simulation studies and two real data analyses.

There are a wide range of psychometric models for analyzing data in educational and psychological surveys. Models including discrete and continuous latent factors have received great attention due to repeated empirical evidence of adequate model fit and success of interpretation that aligns with substantive theory. Among them, a variety of multidimensional item response theory (MIRT) models (Reckase, 2009) have been proposed to account for various multidimensional structures of the latent constructs. Within the family of MIRT models, one of the most studied models is the multidimensional two-parameter logistic model (M2PL) (Reckase, 2009) for dichotomous response, as well as the multidimensional three-parameter logistic model (M3PL) and Multidimensional Four-parameter Logistic Model (M4PL), which are often used when students can guess the answer to the test item correctly (resulting in a lower asymptote in educational measurement) or when the chance of answering an item correctly does not approach 1 (resulting in a higher asymptote in psychopathology measurement). Moreover, the multidimensional graded response model (Cai, 2010) and multidimensional (generalized) partial credit model (Yao & Schwarz, 2006) have also been proposed to handle polytomous items. These models are regarded as extensions of IRT models for dichotomous response in various ways to characterize latent cognitive structures from more complicated datasets.

In the IRT literature, the marginal maximum likelihood estimator (MMLE) is considered to be a consistent and efficient approach for parameter estimation by maximizing the marginal log-likelihood of item parameters (Bock & Aitkin, 1981; Bock et al., 1988). In particular, this estimator is efficient in the asymptotic regime where the test length is small with a rather large number of examinees, which is often the case in applications. However, integrating out the latent ability for marginal likelihood evaluation is notoriously time-consuming, involving multidimensional integrals whose computational complexity grows exponentially with the number of latent dimensions. This results in great difficulties especially when the latent structure is complex or within a domain of large dimension. Several methods have been proposed previously to deal with this computational challenge, such as Gaussian quadrature methods (Bock & Aitkin, 1981; Schilling & Bock, 2005), Laplace approximation methods (Lindstrom & Bates, 1988; Wolfinger & O’connell, 1993; Andersson & Xin, 2021), Metropolis-Hastings Robbins-Monro algorithms (Cai, 2010), stochastic expectation maximization (EM) algorithms (von Davier & Sinharay, 2010; Zhang et al., 2020), and variational approximation methods (Rijmen & Jeon, 2013; Cho et al., 2021; 2022). Among these, Gauss-Hermite quadrature approximation does not scale well to high-dimensional scenarios, and Laplace approximation, which is closely related to variational method (Opper & Archambeau, 2009), may suffer from numerical inaccuracies when dimensions get high or when the likelihood function is in a skewed shape. Additionally, many variants of Laplace approximations, though overcoming some deficiencies, suffer from inflexibility or may be hard to implement (Ormerod & Wand, 2010). The methods based on Monte Carlo simulations such as Metropolis-Hastings Robbins-Monro and stochastic EM algorithms, on the other hand, are more robust but may be computationally inefficient. Notably, in addition to being numerically accurate and computationally efficient, the variational method also provides good interpretability and the variational distributions contain additional useful information (Blei et al., 2017).

Despite the extensive literature on estimation methods for models of dichotomous response (Cai, 2010; Chen et al., 2019; Cho et al., 2021; Feuerstahler & Waller, 2014; Meng et al., 2020), little attention has been given to the efficient and robust estimation of MIRT models for polytomous responses (Bock et al., 1988; Kim & Wilson, 2020). In the paper, we propose a Gaussian variational expectation-maximization algorithm for the multidimensional generalized partial credit model (MGPCM), which is both computationally efficient and numerically stable. The variational method, first proposed in computer science and statistical learning, has since become an efficient approach for large-scale computation in fields such as pattern recognition and document retrieval (Titterington, 2004; Blei & Jordan, 2006). The application of the variational method in analyzing statistical models has also received wide attention (Blei et al., 2017; Ormerod & Wand, 2010). For instance, Ormerod and Wand (2012) adopted a Gaussian variational approximation approach for the estimation of generalized linear mixed effects models. In the field of psychometrics and educational measurement, variational methods have been used for efficient estimation of multidimensional 1PL models (Rijmen & Jeon, 2013) and multidimensional 2/3PL models (Cho et al., 2021), as well as analysis in cognitive diagnostic models (Yamaguchi & Okada, 2020; Oka & Okada 2023). Motivated by Cho et al. (2021), in this paper, we generalize their method to the estimation of MGPCM. This generalization is nontrivial because the MGPCM uses a generalized logit link function that requires a completely new derivation of the variational lower bound to fully enable closed-form parameter update in the EM algorithm.

The rest of the paper is organized as follows. Section 1 provides a brief introduction to the model. In Sect. 2, we introduce the polytomous Gaussian variational expectation-maximization algorithm (pGVEM) and give its derivation. Section 3 evaluates the performance of the pGVEM algorithm compared to the traditional EM algorithm and other existing methods through a comprehensive simulation study. In Sect. 4, we apply the pGVEM algorithm to analyze data from an international educational assessment database and a Big-Five personality assessment. Finally, in Sect. 5, we conclude the paper and suggest potential future research directions.

1. Model Setting

Generalized partial credit model (GPCM) (Muraki, 1992; Embretson & Reise, 2013), also named as a compensatory multidimensional two-parameter partial credit model (M-2PPC) (Yao & Schwarz, 2006), is a popular IRT model for polytomous responses. It allows for the assessment of partial scores for constructed response items and intermediate steps that students have accomplished on the path toward solving the items completely. Suppose we have N examinees and J test items, with the random variable $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} denoting the partial credit of person i’s response to item j. For each item j, $a_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_j$$\end{document} is a discrimination parameter that implies the strength of association between latent trait and item responses. ${\beta }_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{jk}$$\end{document} is a threshold parameter that separates two adjacent response categories. The partial credit model (Masters, 1982) is a special case of GPCM where the discrimination parameter is fixed to be the same across different items. The item response function of GPCM that characterizes the probability of a specific response is given by

(1.1) $\begin{array}{c}Pr(Y_{\mathrm{ij}}=k|{\theta }_i,a_j,{\beta }_{\mathrm{jk}})=\frac{exp\left({\sum }_{r=0}^k,a_j,({\theta }_i-{\beta }_{\mathrm{jr}})\right)}{{\sum }_{v=0}^{K_j}exp\left({\sum }_{r=0}^v,a_j,({\theta }_i-{\beta }_{\mathrm{jr}})\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} \Pr (Y_{ij}=k|\theta _i,a_j,\beta _{jk})=\frac{\exp \left[ \sum _{r=0}^ka_j(\theta _i-\beta _{jr})\right] }{\sum _{v=0}^{K_j}\exp \left[ \sum _{r=0}^va_j(\theta _i-\beta _{jr})\right] }, \end{aligned}$$\end{document}

where $k=0,1,\dots ,K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0,1,\ldots ,K_j-1$$\end{document} and $K_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j$$\end{document} is the number of differential partial credit scores for the $j^{\mathrm{th}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j^{th}$$\end{document} item.

The multidimensional generalized partial credit model (MGPCM) is a natural multidimensional extension of GPCM. The main idea is replacing the uni-dimensional latent ability with a D-dimensional vector, and each dimension represents a facet of the multidimensional construct (i.e., science and literacy). Similarly, the discrimination parameters $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} also become D-dimensional vectors to reflect the discrimination power of item j with respect to each facet (i.e., dimension) of the multidimensional construct $\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} . The threshold parameters stay the same as in uni-dimensional models. Equation (1.1) therefore is updated as follows:

(1.2) $\begin{array}{c}Pr(Y_{\mathrm{ij}}=k|{\theta }_i,a_j,{\beta }_{\mathrm{jk}})=\frac{exp\left\{{\sum }_{r=0}^k(a_j^{\prime }{\theta }_i-{\beta }_{\mathrm{jr}})\right\}}{{\sum }_{v=0}^{K_j-1}exp\left({\sum }_{r=0}^v,(a_j^{\prime }{\theta }_i-{\beta }_{\mathrm{jr}})\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} \Pr (Y_{ij}=k|{\varvec{\theta }}_i,{\varvec{a}}_j,\beta _{jk})=\frac{\exp \{\sum _{r=0}^k({\varvec{a}}_j^\prime {\varvec{\theta }}_i-\beta _{jr})\}}{\sum _{v=0}^{K_j-1}\exp \left\{ \sum _{r=0}^v({\varvec{a}}_j^\prime {\varvec{\theta }}_i-\beta _{jr})\right\} }. \end{aligned}$$\end{document}

Here $a_j^{\prime }{\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{a}}_j^\prime {\varvec{\theta }}_i$$\end{document} indicates the inner product of $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} 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} as $a_j^{\prime }{\theta }_i={\sum }_{d=1}^D{a}_{\mathrm{jd}}{\theta }_{\mathrm{id}}$ \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^\prime {\varvec{\theta }}_i=\sum _{d=1}^Da_{jd}\theta _{id}$$\end{document} where $a_{\mathrm{jd}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jd}$$\end{document} and ${\theta }_{\mathrm{id}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{id}$$\end{document} are the dth component of $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} 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} , respectively. With a slight re-parameterization, we have the following item response function for MGPCM which we will use throughout the paper:

(1.3) $\begin{array}{c}Pr(Y_{\mathrm{ij}}=k|{\theta }_i,a_j,b_{\mathrm{jk}})=\frac{exp(k{a}_j^{\prime }{\theta }_i-b_{\mathrm{jk}})}{{\sum }_{v=0}^{K_j-1}exp(v{a}_j^{\prime }{\theta }_i-b_{\mathrm{jv}})}.\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} \Pr (Y_{ij}=k|{\varvec{\theta }}_i,{\varvec{a}}_j,b_{jk})=\frac{\exp (k{\varvec{a}}_j^{\prime }{\varvec{\theta }}_i-b_{jk})}{\sum _{v=0}^{K_j-1}\exp (v{\varvec{a}}_j^{\prime }{\varvec{\theta }}_i-b_{jv})}. \end{aligned}$$\end{document}

In Equation (1.3), $b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} replaces ${\sum }_{r=0}^k{\beta }_{\mathrm{jr}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{r=0}^k\beta _{jr}$$\end{document} in Equation (1.2) for each $k=0,1,\dots ,K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0,1,\ldots ,K_j-1$$\end{document} . Note that for model identification, we can only have $K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_j-1$$\end{document} estimable threshold parameters, and hence we fix $b_{j0}=0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{j0}=0$$\end{document} .

2. Gaussian Variational Approximation

2.1. Derivation of Algorithm

In this section, we describe the derivation of the GVEM algorithm. In the following, we denote the collection of item parameters by $M_p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_p$$\end{document} being the total number of parameters to be estimated, i.e., $M_p=\{a_j\in R^D,b_{\mathrm{jk}}\in R:j=1,\dots ,J,k=1,\dots ,K_j-1\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_p=\{{\varvec{a}}_j\in {\mathbb {R}}^D,b_{jk}\in {\mathbb {R}}:j=1,\ldots ,J,k=1,\ldots ,K_j-1\}$$\end{document} for MGPCM. As we discussed above, the parameter $b_{j0}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{j0}$$\end{document} is fixed as 0 for all j. To be consistent with the common convention, we assume the latent vector $\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} follows a multivariate normal distribution of $0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{0}}$$\end{document} mean 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}$${\varvec{\Sigma }}_{\theta }$$\end{document} with density function denoted by $\varphi (·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi (\cdot )$$\end{document} . The marginal probability of response vector $Y_i={(Y_{i1},\dots ,Y_{\mathrm{iJ}})}^{\prime }$ \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=(Y_{i1},\ldots ,Y_{iJ})^\prime $$\end{document} for person i is defined as follows

$\begin{array}{cc}Pr\left(Y_i\right|M_p)& ={\int }_{{\theta }_i}\prod _{j=1}^J\prod _{k=0}^{K_j-1}{I}_{(Y_{\mathrm{ij}}=k)}Pr(Y_{\mathrm{ij}}=k|{\theta }_i,M_p)\varphi \left({\theta }_i\right)d{\theta }_i\\ & ={\int }_{{\theta }_i}\prod _{j=1}^J\prod _{k=0}^{K_j-1}{I}_{(Y_{\mathrm{ij}}=k)}\frac{exp(k{a}_j^{\prime }{\theta }_i-b_{\mathrm{jk}})}{{\sum }_{v=0}^{K_j-1}exp(v{a}_j^{\prime }{\theta }_i-b_{\mathrm{jv}})}\varphi \left(\theta \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} \Pr ({\varvec{Y}}_{i}|M_p)&=\int _{{\varvec{\theta }}_i}\prod _{j=1}^J\prod _{k=0}^{K_j-1}I_{(Y_{ij}=k)}\Pr (Y_{ij}=k|{\varvec{\theta }}_i,M_p)\phi ({\varvec{\theta }}_i)d{\varvec{\theta }}_i\\&=\int _{{\varvec{\theta }}_i}\prod _{j=1}^J\prod _{k=0}^{K_j-1}I_{(Y_{ij}=k)}\frac{\exp (k{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jk})}{\sum _{v=0}^{K_j-1}\exp (v{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jv})}\phi ({\varvec{\theta }})d{\varvec{\theta }}_i, \end{aligned}$$\end{document}

where $I_{(Y_{\mathrm{ij}}=k)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{(Y_{ij}=k)}$$\end{document} is an indicator function equal to 1 if $Y_{\mathrm{ij}}=k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij}=k$$\end{document} and zero otherwise, and therefore the marginal log-likelihood function for all responses from the examinees is given by

(2.1) $\begin{array}{c}l\left(M_p\right|Y)=\sum _{i=1}^{N}logPr\left(Y_i\right|M_p)=\sum _{i=1}^{N}log{\int }_{{\theta }_i}\prod _{j=1}^{J}Pr\left(Y_{\mathrm{ij}}\right|{\theta }_i,M_p)\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(M_p|{\varvec{Y}})=\sum _{i=1}^N\log \Pr ({\varvec{Y}}_i|M_p)=\sum _{i=1}^N\log \int _{{\varvec{\theta }}_i}\prod _{j=1}^J\Pr (Y_{ij}|{\varvec{\theta }}_i,M_p)\phi ({\varvec{\theta }}_i)d{\varvec{\theta }}_i, \end{aligned}$$\end{document}

where $Y={(Y_1,\dots ,Y_N)}^{\prime }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Y}}=({\varvec{Y}}_1,\ldots ,{\varvec{Y}}_{N})^\prime $$\end{document} is the $N\times J$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times J$$\end{document} matrix of realized categorical responses.

Following the variational estimation literature (Blei et al., 2017; Cho et al., 2021), we first derive a variational lower bound for MGPCM. We define $KL\left\{p\right(·)\Vert q(·\left)\right\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$KL\{p(\cdot )\Vert q(\cdot )\}$$\end{document} as the Kullback-Leibler divergence of probability distribution p and q. For an arbitrary probability density function $q(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(\cdot )$$\end{document} , the marginal log-likelihood in Equation (2.1) has the following lower bound (Blei et al., 2017):

(2.2) $\begin{array}{cc}& l\left(M_p\right|Y)=\sum _{i=1}^N{\int }_{{\theta }_i}{q}_i\left({\theta }_i\right)d{\theta }_{i}logPr\left(Y_i\right|M_p)\\ & =\sum _{i=1}^N{\int }_{{\theta }_i}log[\frac{P(Y_i,{\theta }_i|M_p)q_i\left({\theta }_i\right)}{P\left({\theta }_i\right|Y_i,M_p)q_i\left({\theta }_i\right)}]q_i\left({\theta }_i\right)d{\theta }_i\\ & =\sum _{i=1}^N[{\int }_{{\theta }_i}[logP(Y_i,{\theta }_i|M_p)]q_i\left({\theta }_i\right)d{\theta }_i\\ & -{\int }_{{\theta }_i}[log{q}_i\left({\theta }_i\right)]q_i\left({\theta }_i\right)d{\theta }_i+KL\{q_i\left({\theta }_i\right)\Vert P\left({\theta }_i\right|Y_i,M_p)\}]\\ & \ge \sum _{i=1}^N{\int }_{{\theta }_i}[logP(Y_i,{\theta }_i|M_p)]q_i\left({\theta }_i\right)d{\theta }_i-\sum _{i=1}^N{\int }_{{\theta }_i}[log{q}_i\left({\theta }_i\right)]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(M_p|{\varvec{Y}})=\sum _{i=1}^N\int _{{\varvec{\theta }}_i}q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i\log \Pr ({\varvec{Y}}_i|M_p)\nonumber \\&\quad =\sum _{i=1}^N\int _{{\varvec{\theta }}_i}\log \Big [\frac{P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)q_i({\varvec{\theta }}_i)}{P({\varvec{\theta }}_i|{\varvec{Y}}_i,M_p)q_i({\varvec{\theta }}_i)}\Big ]q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i\nonumber \\&\quad =\sum _{i=1}^N\bigg {[}\int _{{\varvec{\theta }}_i}\big [\log P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)\big ] q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i\nonumber \\&\quad \quad -\int _{{\varvec{\theta }}_i}\big [\log q_i({\varvec{\theta }}_i)\big ] q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i+KL\Big \{q_i({\varvec{\theta }}_i)\big \Vert P({\varvec{\theta }}_i|{\varvec{Y}}_i,M_p)\Big \}\bigg {]}\nonumber \\&\quad \ge \sum _{i=1}^N\int _{{\varvec{\theta }}_i}\big [\log P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)\big ] q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i-\sum _{i=1}^N\int _{{\varvec{\theta }}_i}\big [\log q_i({\varvec{\theta }}_i)\big ] q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i. \end{aligned}$$\end{document}

The last inequality holds if and only if the KL divergence between the variational distribution $q_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(\cdot )$$\end{document} and the posterior distribution $P(·|Y_i,M_p)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(\cdot |{\varvec{Y}}_i,M_p)$$\end{document} is 0, which indicates $q_i\left({\theta }_i\right)=P\left({\theta }_i\right|Y_i,M_p)$ \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|{\varvec{Y}}_i,M_p)$$\end{document} . In the literature of variational inference, the right-hand side of Equation (2.2) is defined to be the evidence lower bound (Blei et al., 2017), which is equivalent, up to a constant with respect to $q(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(\cdot )$$\end{document} , to the KL divergence between the assumed variational distribution $q(·)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(\cdot )$$\end{document} and the conditional density of the latent variables given the observations.

In the following, we construct an approximation for the marginal maximum likelihood estimator from the evidence lower bound. The primary objective is to identify a suitable distribution that can approximate the posterior distribution $P\left({\theta }_i\right|Y_i,M_p)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{\theta }}_i|{\varvec{Y}}_i,M_p)$$\end{document} . Motivated by this insight, we propose to construct an EM-type algorithm to compute the marginal maximum likelihood estimator. In the E-step, we evaluate the expectation of the complete data log-likelihood, and the expectation is taken with respect to the latent variables ${\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} under its variational probability density function $q_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(\cdot )$$\end{document} :

$\begin{array}{c}\sum _{i=1}^N{\int }_{{\theta }_i}logP(Y_i,{\theta }_i|M_p)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} \sum _{i=1}^N\int _{{\varvec{\theta }}_i}\log P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i, \end{aligned}$$\end{document}

Here the density function $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 chosen to minimize the KL divergence $KL\{q_i({\theta }_i\Vert P({\theta }_i|Y_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\Vert P({\varvec{\theta }}_i|{\varvec{Y}}_i,$$\end{document} $M_p\left)\right\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_p)\}$$\end{document} as the best approximation to the posterior distribution. The second term in the evidence lower bound is left out since it is irrelevant to item parameters. However, a problem with respect to minimizing the KL divergence is that it is hard to find an explicit formula for the posterior 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} with respect to the previous estimated item parameters ${\hat{M}}_p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{M}}_p$$\end{document} , as it involves computing D-dimensional integrals. Numerical methods, such as the Gauss-Hermite approximation, Monte Carlo expectation-maximization, and stochastic expectation-maximization, are often used to provide fast approximation. Herein we adopt the Gaussian variational inference method. It is widely accepted that the posterior distribution of the latent ability $P\left({\theta }_i\right|Y_i,M_p)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{\theta }}_i|Y_i,M_p)$$\end{document} can be approximated by a Gaussian distribution (Chang & Stout, 1993; Wang 2015), and hence we aim to find an optimal $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} in the family of Gaussian distribution while minimizing the KL divergence between $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} and $P\left({\theta }_i\right|Y_i,M_p)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\varvec{\theta }}_i|{\varvec{Y}}_i,M_p)$$\end{document} . Since the posterior distribution can be expressed as

$\begin{array}{c}P\left({\theta }_i\right|Y_i,M_p)=\frac{P(Y_i,{\theta }_i|M_p)}{P\left(Y_i\right|M_p)},\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({\varvec{\theta }}_i|{\varvec{Y}}_i, M_p)=\frac{P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)}{P({\varvec{Y}}_i|M_p)}, \end{aligned}$$\end{document}

we only need to evaluate $P(Y_i,{\theta }_i|M_p)$ \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,{\varvec{\theta }}_i|M_p)$$\end{document} to find a proper $q_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(\cdot )$$\end{document} as

$\begin{array}{c}KL\{q_i\left({\theta }_i\right)\left|\right|P\left({\theta }_i\right|Y_i,M_p)\}=KL\{q_i\left({\theta }_i\right)\left|\right|P({\theta }_i,Y_i|M_p)\}+C.\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)||P({\varvec{\theta }}_i|{\varvec{Y}}_i,M_p)\}=KL\{q_i({\varvec{\theta }}_i)||P({\varvec{\theta }}_i,{\varvec{Y}}_i|M_p)\}+C. \end{aligned}$$\end{document}

Under the setting of MGPCM, the logarithm of joint distribution function 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

(2.3) $\begin{array}{cc}& logP(Y_i,{\theta }_i|M_p)=logP\left(Y_i\right|{\theta }_i,M_p)+log\varphi \left({\theta }_i\right)\\ & =\sum _{j=1}^J\left(\sum _{k=0}^{K_j-1},I_{(Y_{\mathrm{ij}}=k)},\left(k,a_j^{\prime },{\theta }_i,-,b_{\mathrm{jk}},-,log,\left(\sum _{v=0}^{K_j-1},exp,(v{a}_j^{\prime }{\theta }_i-b_{\mathrm{jv}})\right)\right)\right)+log\varphi \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({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)=\log P({\varvec{Y}}_i|{\varvec{\theta }}_i,M_p)+\log \phi ({\varvec{\theta }}_i)\nonumber \\&\quad =\sum _{j=1}^J\left\{ \sum _{k=0}^{K_j-1}I_{(Y_{ij}=k)}\left[ k{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jk}-\log \left( \sum _{v=0}^{K_j-1}\exp (v{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jv})\right) \right] \right\} +\log \phi ({\varvec{\theta }}_i). \end{aligned}$$\end{document}

The nonlinear softmax function, defined by $f_v\left(x\right)=exp\left(x_v\right)/[{\sum }_{k=1}^{n}exp\left(x_k\right)]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_v({\varvec{x}})=\exp (x_v)/[\sum _{k=1}^n\exp (x_k)]$$\end{document} for an n-dimensional vector $x$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{x}}$$\end{document} , is the main cause of the intractability of integral. To overcome this problem, a variational lower bound based on the approximation to the softmax function is proposed and by augmenting Equation (2.3) with variational parameters, the evidence lower bound can be computed explicitly without resorting to numeric integration. Among the many approximations for the softmax function, we adopt a One-Versus-Each bound (Tisais, 2016), which well approximates the softmax function and captures the model features. We start with the following inequality:

(2.4) $\begin{array}{c}{f}_v\left(x\right)=\frac{e^{x_v}}{{\sum }_{k=1}^n{e}^{x_k}}\ge \prod _{k=1,k\ne v}^n\frac{e^{x_v}}{e^{x_v}+e^{x_k}}=2\prod _{k=1}^n\frac{e^{x_v}}{e^{x_v}+e^{x_k}}.\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} f_v({\varvec{x}})=\frac{e^{x_v}}{\sum _{k=1}^ne^{x_k}}\ge \prod _{k=1,k\ne v}^n\frac{e^{x_v}}{e^{x_v}+e^{x_k}}=2\prod _{k=1}^n\frac{e^{x_v}}{e^{x_v}+e^{x_k}}. \end{aligned}$$\end{document}

Denote by $k_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{ij}$$\end{document} the realized partial credit score for the response of the ith person for the jth item. Then by applying bound (2.4) to (2.3), we have

$\begin{array}{cc}logP(Y_i,{\theta }_i|M_p)& =log\varphi \left({\theta }_i\right)+\sum _{j=1}^J\left(\sum _{k=0}^{K_j-1},I_{(Y_{\mathrm{ij}}=k)},log,\frac{exp\left(x_{\mathrm{ijk}}\right)}{{\sum }_{v=0}^{K_j-1}exp\left(x_{\mathrm{ijv}}\right)}\right)\\ & \ge log\varphi \left({\theta }_i\right)+\sum _{j=1}^J\left(\sum _{k=0}^{K_j-1},I_{(Y_{\mathrm{ij}}=k)},\left(log,2,+,\sum _{v=0}^{K_j-1},log,\frac{exp\left(x_{\mathrm{ijk}}\right)}{exp\left(x_{\mathrm{ijv}}\right)+exp\left(x_{\mathrm{ijk}}\right)}\right)\right)\\ & =log\varphi \left({\theta }_i\right)+\sum _{j=1}^J\left(log,2,-,\sum _{k=0}^{K_j-1},log,\left(1,+,exp,(,x_{\mathrm{ijk}},-,x_{ij{k}_{\mathrm{ij}}},)\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} \log P({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)&=\log \phi ({\varvec{\theta }}_i)+\sum _{j=1}^J\left[ \sum _{k=0}^{K_j-1}I_{(Y_{ij}=k)}\log \frac{\exp (x_{ijk})}{\sum _{v=0}^{K_j-1}\exp (x_{ijv})}\right] \\&\ge \log \phi ({\varvec{\theta }}_i){+}\sum _{j=1}^J \left\{ \sum _{k=0}^{K_j-1}I_{(Y_{ij}=k)}\left[ \log 2{+}\sum _{v=0}^{K_j-1}\log \frac{\exp (x_{ijk})}{\exp (x_{ijv})+\exp (x_{ijk})}\right] \right\} \\&=\log \phi ({\varvec{\theta }}_i)+\sum _{j=1}^J\left\{ \log 2-\sum _{k=0}^{K_j-1}\log \left[ 1+\exp (x_{ijk}-x_{ijk_{ij}})\right] \right\} . \end{aligned}$$\end{document}

Here we denote $x_{\mathrm{ijk}}=k{a}_j^{\prime }{\theta }_i-b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ijk}=k{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jk}$$\end{document} for short. We wish to draw attention to our selection of the “One-Versus-Each bound.” It can be established for (2.4) that a strict inequivalence holds true in all cases, except for the exceptional circumstance $x_v/x_k\to \infty$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x_v}/{x_k}\rightarrow \infty $$\end{document} for all $k\ne v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k\ne v$$\end{document} with at most one exception. Additionally, the approximation is closest when $x_v$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_v$$\end{document} is among the largest of all $x_k$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_k$$\end{document} . The idea of maximum likelihood estimation indicates that, when partial credit score $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} is recorded as $k_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{ij}$$\end{document} , $k_{\mathrm{ij}}{a}_j^{\prime }{\theta }_i-b_{j{k}_{\mathrm{ij}}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{ij}{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jk_{ij}}$$\end{document} is the most likely to be the largest among all $v{a}_j^{\prime }{\theta }_i-b_{\mathrm{jv}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jv}$$\end{document} for $v=0,1,\dots ,K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v=0,1,\ldots ,K_j-1$$\end{document} . Therefore the feature of the One-Versus-Each bound does fit well as an approximation to the marginal maximum likelihood estimator.

Logistic sigmoid function (2.4) can be further approximated by a local variational approach:

$\begin{array}{cc}logP(Y_i,{\theta }_i|M_p)& \ge log\varphi \left({\theta }_i\right)+\sum _{j=1}^J\left(log,2,-,\sum _{k=0}^{K_j-1},\eta ,\left({\xi }_{\mathrm{ijk}}\right),\left({(x_{\mathrm{ijk}}-x_{ij{k}_{\mathrm{ij}}})}^2,-,{{\xi }_{\mathrm{ijk}}}^2\right)\right)\\ & \left(-,\sum _{k=0}^{K_j-1},\frac{1}{2},(x_{\mathrm{ijk}}-x_{ij{k}_{\mathrm{ij}}}-{\xi }_{\mathrm{ijk}}),-,\sum _{k=0}^{K_j-1},log,(1+e^{{\xi }_{\mathrm{ijk}}})\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({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)&\ge \log \phi ({\varvec{\theta }}_i)+\sum _{j=1}^J\left\{ \log 2-\sum _{k=0}^{K_j-1}\eta (\xi _{ijk}) \left[ (x_{ijk}-x_{ijk_{ij}})^2-{\xi _{ijk}}^2\right] \right. \\&\quad \left. -\sum _{k=0}^{K_j-1}\frac{1}{2}(x_{ijk}-x_{ijk_{ij}}-\xi _{ijk})- \sum _{k=0}^{K_j-1}\log (1+e^{\xi _{ijk}}) \right\} , \end{aligned}$$\end{document}

where $\xi ={\left\{{\xi }_{\mathrm{ijk}}\right\}}_{i,j,k}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}=\{\xi _{ijk}\}_{i,j,k}$$\end{document} are called variational parameters, which will be iteratively updated together with item parameters in the M-step. Here the function $\eta \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}$$\eta (x)$$\end{document} is defined as $(e^x-1)/\left[4x(e^x+1)\right]$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${(e^x-1)}/{[4x(e^x+1)]}$$\end{document} (Jaakkola & Jordan, 2000). We use this local variational approximation for a suitable expectation of $logP(Y_i,{\theta }_i|M_p)$ \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{Y}}_i,{\varvec{\theta }}_i|M_p)$$\end{document} (i.e., given below in Equation (3.5)) that can be written as a quadratic form with respect to ${\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} . This will facilitate the selection of $q_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(\cdot )$$\end{document} in the family of Gaussian distributions.

Next we substitute $x_{\mathrm{ijk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{ijk}$$\end{document} by $k{a}_j^{\prime }{\theta }_i-b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k{\varvec{a}}_j^\prime {\varvec{\theta }}_i-b_{jk}$$\end{document} and write the above lower bound of joint distribution function as

$\begin{array}{cc}logP(Y_i,{\theta }_i|M_p)& \ge \sum _{j=1}^J\left(-,,\sum _{k=0}^{K_j-1},\left(\eta ,\left({\xi }_{\mathrm{ijk}}\right),{(k-k_{\mathrm{ij}})}^2,{\theta }_i^{\prime },a_j,a_j^{\prime },{\theta }_i,-,2,(k-k_{\mathrm{ij}}),\eta ,\left({\xi }_{\mathrm{ijk}}\right),(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}}),a_j^{\prime },{\theta }_i\right)\right)\\ & \left(\left(+,\frac{1}{2},(k-k_{\mathrm{ij}}),a_j^{\prime },{\theta }_i,+,\eta ,\left({\xi }_{\mathrm{ijk}}\right),{(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})}^2,-,\frac{1}{2},(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})\right)\right)\\ & \left(\left(-,\eta ,\left({\xi }_{\mathrm{ijk}}\right),{{\xi }_{\mathrm{ijk}}}^2,-,\frac{1}{2},{\xi }_{\mathrm{ijk}},+,log,(1+e^{{\xi }_{\mathrm{ijk}}})\right),+,log,2\right)+log\varphi \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({\varvec{Y}}_i,{\varvec{\theta }}_i|M_p)&\ge \sum _{j=1}^J\left\{ {-}\!\sum _{k=0}^{K_j-1} \left[ \eta (\xi _{ijk})(k{-}k_{ij})^2{\varvec{\theta }}_i^\prime {\varvec{a}}_j{\varvec{a}}_j^\prime {\varvec{\theta }}_i{-}2(k-k_{ij})\eta (\xi _{ijk}) (b_{jk}{-}b_{jk_{ij}}){\varvec{a}}_j^\prime {\varvec{\theta }}_i \right. \right. \\&\quad \left. \left. + \frac{1}{2}(k-k_{ij}){\varvec{a}}_j^\prime {\varvec{\theta }}_i+\eta (\xi _{ijk}) (b_{jk}-b_{jk_{ij}})^2-\frac{1}{2}(b_{jk}-b_{jk_{ij}})\right. \right. \\&\quad \left. \left. -\eta (\xi _{ijk}){\xi _{ijk}}^2-\frac{1}{2}\xi _{ijk}+\log (1+e^{\xi _{ijk}})\right] +\log 2\right\} +\log \phi ({\varvec{\theta }}_i), \end{aligned}$$\end{document}

and therefore the expectation of the log-likelihood, which needs computing in the E-step, takes the following form:

(2.5) $\begin{array}{cc}E(M_p,\xi )& :={\int }_{{\theta }_i}logP\left(Y_i\right|{\theta }_i,M_p)q_i\left({\theta }_i\right)d{\theta }_i+{\int }_{{\theta }_i}log\varphi \left({\theta }_i\right)q_i\left({\theta }_i\right)d{\theta }_i\\ & \ge {\int }_{{\theta }_i}\sum _{j=1}^J\{log2-\sum _{k=0}^{K_j-1}[\eta \left({\xi }_{\mathrm{ijk}}\right){(k-k_{\mathrm{ij}})}^2{\theta }_i^{\prime }{a}_j{a}_j^{\prime }{\theta }_i-2(k-k_{\mathrm{ij}})\eta \left({\xi }_{\mathrm{ijk}}\right)(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})a_j^{\prime }{\theta }_i\\ & +\frac{1}{2}(k-k_{\mathrm{ij}})a_j^{\prime }{\theta }_i+\eta \left({\xi }_{\mathrm{ijk}}\right){(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})}^2-\frac{1}{2}(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})\\ & -\eta \left({\xi }_{\mathrm{ijk}}\right){{\xi }_{\mathrm{ijk}}}^2-\frac{1}{2}{\xi }_{\mathrm{ijk}}+log(1+e^{{\xi }_{\mathrm{ijk}}})]\}{q}_i\left({\theta }_i\right)d{\theta }_i+{\int }_{{\theta }_i}log\varphi \left({\theta }_i\right)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} E(M_p,{\varvec{\xi }})&:=\int _{{\varvec{\theta }}_i}\log P({\varvec{Y}}_i|{\varvec{\theta }}_i,M_p)q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i+\int _{{\varvec{\theta }}_i}\log \phi ({\varvec{\theta }}_i)q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i\nonumber \\&\ge \int _{{\varvec{\theta }}_i}\sum _{j=1}^J\Bigg \{\log 2{-}\!\sum _{k=0}^{K_j-1}\Big [\eta (\xi _{ijk})(k{-}k_{ij})^2{\varvec{\theta }}_i^\prime {\varvec{a}}_j{\varvec{a}}_j^\prime {\varvec{\theta }}_i{-}2(k{-}k_{ij})\eta (\xi _{ijk})(b_{jk}{-}b_{jk_{ij}}){\varvec{a}}_j^\prime {\varvec{\theta }}_i\nonumber \\&\quad +\frac{1}{2}(k-k_{ij}){\varvec{a}}_j^\prime {\varvec{\theta }}_i+\eta (\xi _{ijk})(b_{jk}-b_{jk_{ij}})^2-\frac{1}{2}(b_{jk}-b_{jk_{ij}})\nonumber \\&\quad -\eta (\xi _{ijk}){\xi _{ijk}}^2-\frac{1}{2}\xi _{ijk}+\log (1+e^{\xi _{ijk}})\Big ]\Bigg \}q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i+\int _{{\varvec{\theta }}_i}\log \phi ({\varvec{\theta }}_i)q_i({\varvec{\theta }}_i)d{\varvec{\theta }}_i. \end{aligned}$$\end{document}

For a minimized KL divergence, $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(\theta _i)$$\end{document} is selected as follows:

$\begin{array}{cc}& log{q}_i\left({\theta }_i\right)\propto \sum _{j=1}^J\sum _{k=0}^{K_j-1}\{(k-k_{\mathrm{ij}})[2\eta \left({\xi }_{\mathrm{ijk}}\right)(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})-0.5]a_j^{\prime }{\theta }_i-\eta \left({\xi }_{\mathrm{ijk}}\right){(k-k_{\mathrm{ij}})}^2{\theta }_i^{\prime }{a}_j{a}_j^{\prime }{\theta }_i\}\\ & -\frac{{\theta }_i^{\prime }{\Sigma }_{\theta }^{-1}{\theta }_i}{2}.\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 q_i({\varvec{\theta }}_i)\propto \sum _{j=1}^J\sum _{k=0}^{K_j-1}\Big \{(k{-}k_{ij})\big [2\eta (\xi _{ijk})(b_{jk}{-}b_{jk_{ij}}){-}0.5\big ]{\varvec{a}}_j^\prime {\varvec{\theta }}_i{-}\eta (\xi _{ijk})(k{-}k_{ij})^2{\varvec{\theta }}_i^\prime {\varvec{a}}_j{\varvec{a}}_j^\prime {\varvec{\theta }}_i\Big \}\\ {}{} & {} \quad -\frac{{\varvec{\theta }}_i^\prime {\varvec{\Sigma }}_\theta ^{-1}{\varvec{\theta }}_i}{2}. \end{aligned}$$\end{document}

As the choice of $q_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(\cdot )$$\end{document} has been confined in the Gaussian family, it suffices to give the update for the mean and covariance matrix:

(2.6) $\begin{array}{cc}{\mu }_i=& {\Sigma }_i\times \sum _{j=1}^J\sum _{k=0}^{K_j-1}(k-k_{\mathrm{ij}})[2\eta \left({\xi }_{\mathrm{ijk}}\right)(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})-\frac{1}{2}]a_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} {\varvec{\mu }}_i= & {} {\varvec{\Sigma }}_i\times \sum _{j=1}^J\sum _{k=0}^{K_j-1}(k-k_{ij})\Big [2\eta (\xi _{ijk})(b_{jk}-b_{jk_{ij}})-\frac{1}{2}\Big ]{\varvec{a}}_{j}; \end{aligned}$$\end{document}

(2.7) $\begin{array}{cc}{\Sigma }_i^{-1}=& {\Sigma }_{\theta }^{-1}+2\sum _{j=1}^J\sum _{k=0}^{K_j-1}\eta \left({\xi }_{\mathrm{ijk}}\right){(k-k_{\mathrm{ij}})}^2{a}_j{a}_j^{\prime }.\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 }}_i^{-1}= & {} {{\varvec{\Sigma }}}_\theta ^{-1}+2\sum _{j=1}^J\sum _{k=0}^{K_j-1}\eta (\xi _{ijk})(k-k_{ij})^2{\varvec{a}}_j{\varvec{a}}_j^\prime . \end{aligned}$$\end{document}

In each iteration, the item parameters and variational parameters ${\xi }_{\mathrm{ijk}},a_j,b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ijk},{\varvec{a}}_j,b_{jk}$$\end{document} are obtained from the previous M-step and taken as the initial value if it is the first iteration.

For the M-step, the item parameters are chosen to maximize the above expectation of the lower bound obtained by plugging Equation (2.6) and (2.7) into Equation (2.5):

(2.8) $\begin{array}{cc}E(M_p,\xi )& \ge \sum _{i=1}^N\sum _{j=1}^{J}log2+\sum _{i=1}^N\sum _{j=1}^J\sum _{k=0}^{K_j-1}\{-\eta \left({\xi }_{\mathrm{ijk}}\right){(k-k_{\mathrm{ij}})}^2{a}_j^{\prime }[{\Sigma }_i^{\left(t\right)}+\left({\mu }_i^{\left(t\right)}\right){\left({\mu }_i^{\left(t\right)}\right)}^{\prime }]a_j\\ & +(k-k_{\mathrm{ij}})[2\eta \left({\xi }_{\mathrm{ijk}}\right)(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})-\frac{1}{2}]a_j^{\prime }{\mu }_i^{\left(t\right)}-\eta \left({\xi }_{\mathrm{ijk}}\right){(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})}^2+\frac{1}{2}(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})\\ & +\eta \left({\xi }_{\mathrm{ijk}}\right){{\xi }_{\mathrm{ijk}}}^2+\frac{1}{2}{\xi }_{\mathrm{ijk}}-log(1+e^{{\xi }_{\mathrm{ijk}}})\}-\frac{N}{2}log\left|{\Sigma }_{\theta }^{\left(t\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} E(M_p,{\varvec{\xi }})&\ge \sum _{i=1}^N\sum _{j=1}^J\log 2+\sum _{i=1}^N\sum _{j=1}^J\sum _{k=0}^{K_j-1}\Bigg \{ -\eta (\xi _{ijk})(k-k_{ij})^2{\varvec{a}}_{j}^\prime \Big [{\varvec{\Sigma }}_i^{(t)}+({\varvec{\mu }}_i^{(t)})({\varvec{\mu }}_i^{(t)})^\prime \Big ]{\varvec{a}}_{j}\nonumber \\&\quad +(k-k_{ij})[2\eta (\xi _{ijk})(b_{jk}{-}b_{jk_{ij}}){-}\frac{1}{2}]{\varvec{a}}_j^\prime {\varvec{\mu }}_i^{(t)}{-}\eta (\xi _{ijk})(b_{jk}{-}b_{jk_{ij}})^2{+}\frac{1}{2}(b_{jk}{-}b_{jk_{ij}})\nonumber \\&\quad +\eta (\xi _{ijk}){\xi _{ijk}}^2+\frac{1}{2}\xi _{ijk}-\log (1+e^{\xi _{ijk}}) \Bigg \}-\frac{N}{2}\log |{\varvec{\Sigma }}_{\theta }^{(t)}| \end{aligned}$$\end{document}

(2.9) $\begin{array}{cc}& -\sum _{i=1}^N\frac{1}{2}Tr\{{\left({\Sigma }_{\theta }^{\left(t\right)}\right)}^{-1}[{\Sigma }_i^{\left(t\right)}+\left({\mu }_i^{\left(t\right)}\right){\left({\mu }_i^{\left(t\right)}\right)}^{\prime }]\\ & :=\underline{E}(M_p,\xi ).\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 -\sum _{i=1}^N\frac{1}{2}Tr\{({\varvec{\Sigma }}_{\theta }^{(t)})^{-1}[{\varvec{\Sigma }}_i^{(t)}+({\varvec{\mu }}_i^{(t)})({\varvec{\mu }}_i^{(t)})^\prime ]\nonumber \\&{:}{=}\, \underline{E}(M_p,{\varvec{\xi }}). \end{aligned}$$\end{document}

Through maximizing the lower bound $\underline{E}(M_p,\xi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{E}(M_p,{\varvec{\xi }})$$\end{document} , we can derive a new set of item parameters that could potentially maximize the left-hand side. Updating the variational parameters helps to prevent the iteration from leading to a smaller value of the target expectation by shrinking the inequality too much when the right-hand side is maximized. The efficiency of this majorization-maximization approach depends on the goodness of fit of the adopted softmax bound.

To maximize the lower bound on the expectation concerning the item parameters $M_p$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_p$$\end{document} and variational parameters $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\xi }}$$\end{document} , we employ a Gauss-Seidel scheme to handle the nonlinear terms regarding the parameters. Each iterative update uses the most recently updated copies of the parameters. The update is given as follows.

For each $j=1,\dots ,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,\ldots , J$$\end{document} ,

(2.10) $\begin{array}{cc}{a}_j=& \frac{1}{2}{\left(\sum _{i=1}^N,\sum _{k=0}^{K_j-1},\eta ,\left({\xi }_{\mathrm{ijk}}\right),{(k-k_{\mathrm{ij}})}^2,({\Sigma }_i+{\mu }_i{\mu }_i^{\prime })\right)}^{-1}\\ & \left(\sum _{i=1}^N,\sum _{k=0}^{K_j-1},(k-k_{\mathrm{ij}}),\left(2,\eta ,\left({\xi }_{\mathrm{ijk}}\right),(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}}),-,\frac{1}{2}\right),{\mu }_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} {\varvec{a}}_j= & {} \frac{1}{2}\left[ \sum _{i=1}^N\sum _{k=0}^{K_j-1}\eta (\xi _{ijk})(k-k_{ij})^2({\varvec{\Sigma }}_i+{\varvec{\mu }}_i{\varvec{\mu }}_i^\prime )\right] ^{-1}\nonumber \\{} & {} \left\{ \sum _{i=1}^N\sum _{k=0}^{K_j-1}(k-k_{ij})\left[ 2\eta (\xi _{ijk})(b_{jk}-b_{jk_{ij}})-\frac{1}{2}\right] {\varvec{\mu }}_i\right\} , \end{aligned}$$\end{document}

with ${\xi }_{\mathrm{ijk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ijk}$$\end{document} , $b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} from the last iteration or initialization and ${\mu }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\mu }}_i$$\end{document} , ${\Sigma }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Sigma }}_i$$\end{document} from E-step. Next for the threshold parameters, for each $j=1,\dots ,J,k=1,\cdots K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,J,k=1,\cdots K_j-1$$\end{document} ,

(2.11) $\begin{array}{c}{b}_{\mathrm{jk}}=\frac{{\sum }_{i=1}^N[{\text{B}}_1(i,j,k)I_{(k\ne k_{\mathrm{ij}})}+I_{(k=k_{\mathrm{ij}})}{\sum }_{v=0,v\ne k}^{K_j-1}{\text{B}}_2(i,j,v,k)]}{2{\sum }_{i=1}^N(\eta \left({\xi }_{\mathrm{ijk}}\right)I_{(k\ne k_{\mathrm{ij}})}+I_{(k=k_{\mathrm{ij}})}{\sum }_{v=0,v\ne k}^{K_j-1}\eta \left({\xi }_{\mathrm{ijv}}\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} b_{jk}=\frac{\sum _{i=1}^N[\text{ B}_1(i,j,k)I_{(k\ne k_{ij})}+I_{(k=k_{ij})}\sum _{v=0,v\ne k}^{K_j-1}\text{ B}_2(i,j,v,k)]}{2\sum _{i=1}^N(\eta (\xi _{ijk})I_{(k\ne k_{ij})}+I_{(k=k_{ij})}\sum _{v=0,v\ne k}^{K_j-1}\eta (\xi _{ijv}))}, \end{aligned}$$\end{document}

where

$\begin{array}{cc}{\text{B}}_1(i,j,k)& =2\eta \left({\xi }_{\mathrm{ijk}}\right)(k-k_{\mathrm{ij}})a_j^{\prime }{\mu }_i+0.5+2\eta \left({\xi }_{\mathrm{ijk}}\right)b_{j{k}_{\mathrm{ij}}};\\ {\text{B}}_2(i,j,v,k)& =-2\eta \left({\xi }_{\mathrm{ijk}}\right)(v-k)a_j^{\prime }{\mu }_i-0.5+2\eta \left({\xi }_{\mathrm{ijv}}\right)b_{\mathrm{jv}}.\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{ B}_1(i,j,k)&=2\eta (\xi _{ijk})(k-k_{ij}){\varvec{a}}_j^\prime {\varvec{\mu }}_i+0.5+2\eta (\xi _{ijk})b_{jk_{ij}};\\\text{ B}_2(i,j,v,k)&=-2\eta (\xi _{ijk})(v-k){\varvec{a}}_j^\prime {\varvec{\mu }}_i-0.5+2\eta (\xi _{ijv})b_{jv}. \end{aligned}$$\end{document}

Here $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} are from the previous step and ${\mu }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\mu }}_i$$\end{document} , ${\Sigma }_i$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\Sigma }}_i$$\end{document} are from the previous E-step. Finally for the variational parameters ${\xi }_{\mathrm{ijk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{ijk}$$\end{document} , for each $i=1,\dots ,I,j=1,\dots ,J,k=0,\cdots K_j-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,\ldots ,I,j=1,\ldots , J,k=0,\cdots K_j-1$$\end{document} , we have update as

(2.12) $\begin{array}{c}{\xi }_{\mathrm{ijk}}^2={[(k-k_{\mathrm{ij}})a_j^{\prime }{\mu }_i-(b_{\mathrm{jk}}-b_{j{k}_{\mathrm{ij}}})]}^2+{(k-k_{\mathrm{ij}})}^2{a}_j^{\prime }{\Sigma }_i{a}_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} \xi _{ijk}^2=[(k-k_{ij}){\varvec{a}}_j^\prime {\varvec{\mu }}_i-(b_{jk}-b_{jk_{ij}})]^2+(k-k_{ij})^2{\varvec{a}}_j^\prime {\varvec{\Sigma }}_i{\varvec{a}}_j. \end{aligned}$$\end{document}

with all other parameters obtained from the latest updates.

In the exploratory analysis where we do not have any prior information on the item factor loadings, so the assumed 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 _\theta $$\end{document} is fixed as $I_D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_D$$\end{document} and later proper rotations (Browne, 2001) are imposed to allow the factors to be correlated and thus allow for analysis of latent structures. But for confirmatory factor analysis, we update the covariance as

(2.13) $\begin{array}{c}{\Sigma }_{\theta }=\frac{1}{N}\sum _{i=1}^N({\Sigma }_i+{\mu }_i{\mu }_i^{\prime }).\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 }}_\theta =\frac{1}{N}\sum _{i=1}^N({\varvec{\Sigma }}_i+{\varvec{\mu }}_i{\varvec{\mu }}_i^\prime ). \end{aligned}$$\end{document}

and scale its diagonal entries to be 1.

2.2. Standard Error Estimation

Computing standard errors (SEs) of the item parameter estimates is crucial for various applications, such as multidimensional computerized adaptive testing, item parameter calibration as well as differential item functioning. Challenges arise in estimating SEs when dealing with a high-dimensional latent domain and polytomous responses as in MGPCM. The commonly used method for estimating SEs is based on the approximated Fisher’s information matrix. However, taking the inverse of a prohibitively large information matrix (due to high dimensions and long test length) can be unstable when the sample size of the examinees is not large enough. An alternative numerical approximation using Gaussian quadrature in EM estimation has been proposed (Cagnone & Monari, 2013), but it is computationally expensive and sensitive to dimensionality. The supplemented expectation-maximization (SEM) algorithm has also been developed in the IRT literature (e.g., Tian et al., 2013). However, in pilot simulations we found that none of these methods are capable of providing stable estimations, especially when the dimension D and the number of categories K are large. Therefore, to estimate the standard errors of item parameters under the pGVEM framework for MGPCM, we adopt a bootstrap approach that uses a resampling procedure. Bootstrap is an efficient alternative when the standard SEs estimation is mathematically intractable (Efron & Tibshirani, 1986). The resampling procedure avoids the direct computation of SEs.

The bootstrap procedure in the pGVEM framework is implemented as follows. First we simulated B bootstrap datasets based on ${\hat{M}}_p={\{{\hat{a}}_j,{\hat{b}}_j\}}_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{M}}_p= \{\hat{{\varvec{a}}}_{j},\hat{{\varvec{b}}}_{j}\}_j$$\end{document} estimated from the pGVEM scheme. Then we apply the pGVEM method to estimate the item parameters for each of the bootstrap datasets, denoted by ${\hat{M}}_p^{\left(1\right)},\dots ,{\hat{M}}_p^{\left(B\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{M}}_p^{(1)},\ldots ,{\hat{M}}_p^{(B)}$$\end{document} . The standard errors are estimated by

$\begin{array}{c}{\widehat{\mathrm{SE}}}_v=\sqrt{\frac{1}{B-1}\sum _{i=1}^B{({\hat{v}}^{\left(i\right)}-\hat{v})}^2},\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} {\widehat{SE}}_{v}=\sqrt{\frac{1}{B-1}\sum _{i=1}^B({\hat{v}}^{(i)}-{\hat{v}})^2}, \end{aligned}$$\end{document}

where v denotes item parameter $a_{\mathrm{jr}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jr}$$\end{document} or $b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} and ${\hat{v}}^{\left(i\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{v}}^{(i)}$$\end{document} is its ith bootstrap estimate. Given that our objective is to estimate SEs rather than the distributions of the estimators, in our study, we take the number of bootstrap samples to be 50, which generates stable results numerically.

2.3. Determining Latent Dimension

In this section, we discuss how to select the appropriate number of latent dimensions. We propose to use the information criterion such as AIC or BIC to compare the model fit with different dimensions. In the MGPCM, direct computation of the residual sum of squares is costly. So we adopt the modified version of the information criterion where the expectation is replaced by its lower bound given by (2.9):

(2.14) $\begin{array}{cc}AI{C}^{\ast }& =2(\Vert \hat{A}{\Vert }_0+\Vert \hat{B}{\Vert }_0+\Vert \text{nondiag}\left({\hat{\Sigma }}_{\theta }\right){\Vert }_0/2)-2\underline{E}({\hat{M}}_p,\hat{\xi }),\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} AIC^{*}&=2(\Vert \hat{{\varvec{A}}}\Vert _0+\Vert \hat{{\varvec{B}}}\Vert _0+\Vert \text{ nondiag }(\hat{{\varvec{\Sigma }}}_\theta )\Vert _0/2)-2\underline{E}({\hat{M}}_p,\hat{{\varvec{\xi }}}), \end{aligned}$$\end{document}

(2.15) $\begin{array}{cc}BI{C}^{\ast }& =log\left(N\right)(\Vert \hat{A}{\Vert }_0+\Vert \hat{B}{\Vert }_0+\Vert \text{nondiag}\left({\hat{\Sigma }}_{\theta }\right){\Vert }_0/2)-2\underline{E}({\hat{M}}_p,\hat{\xi }),\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} BIC^{*}&=\log (N)(\Vert \hat{{\varvec{A}}}\Vert _0+\Vert \hat{{\varvec{B}}}\Vert _0+\Vert \text{ nondiag }(\hat{{\varvec{\Sigma }}}_\theta )\Vert _0/2)-2\underline{E}({\hat{M}}_p,\hat{{\varvec{\xi }}}), \end{aligned}$$\end{document}

where $\hat{A}=({\hat{a}}_1,\dots ,{\hat{a}}_J)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{A}}}=(\hat{{\varvec{a}}}_1,\ldots ,\hat{{\varvec{a}}}_J)$$\end{document} and $\hat{B}=({\hat{b}}_1,\dots ,{\hat{b}}_J)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{B}}}=(\hat{{\varvec{b}}}_1,\ldots ,\hat{{\varvec{b}}}_J)$$\end{document} are assembled matrices of the discrimination and threshold parameters, respectively, $\text{nondiag}\left({\hat{\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}$$\text{ nondiag }(\hat{{\varvec{\Sigma }}}_\theta )$$\end{document} denotes the nondiagonal entries of ${\hat{\Sigma }}_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\Sigma }}}_\theta $$\end{document} , and the zero norm $\left|\right|·{\left|\right|}_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$||\cdot ||_0$$\end{document} counts the number of nonzero entries of the assembled matrix. Here note that the term $\Vert \hat{B}{\Vert }_0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \hat{{\varvec{B}}}\Vert _0$$\end{document} does not increase with dimension D and it denotes the number of all effective b parameters. In addition, since the covariance matrix ${\hat{\Sigma }}_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\Sigma }}}_{\theta }$$\end{document} is symmetric with unit diagonal entries, we count the effective number of parameters in ${\hat{\Sigma }}_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{{\varvec{\Sigma }}}_{\theta }$$\end{document} as $\Vert \text{nondiag}\left({\hat{\Sigma }}_{\theta }\right){\Vert }_0/2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \text{ nondiag }(\hat{{\varvec{\Sigma }}}_\theta )\Vert _0/2$$\end{document} . The major advantage of the proposed criteria is that the lower bound of expectation is readily obtained with the updated item parameters and variation parameters, with no extra computation cost.

3. Simulation Studies

3.1. Study I

We conducted simulation studies to compare the empirical performance of the proposed pGVEM algorithm with the EM algorithm with fixed quadrature, Metropolis-Hastings Robbins-Monro (MHRM) algorithm (Cai, 2010), and the stochastic EM (StEM) for MGPCM, which are implemented in the R package ‘mirt’ (Chalmers, 2012), in terms of mean squared error and bias of the estimation together with their computational time. The simulations were conducted in the exploratory factor analysis (EFA) scenario, where no constraints on the item factor loading structure were imposed during the analysis. EFA is generally more computationally challenging than confirmatory factor analysis. In our study, we assumed that the latent covariance matrix was $I_D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_D$$\end{document} to remove scale and rotational indeterminacy, and no further assumptions on the structure of the loading matrix A were made during the analysis.

After estimating the parameters by our pGVEM algorithm, we performed proper oblique rotation to allow the factors to be correlated. Many methods, including varimax, direct oblimin, quartimax, equamax, and promax (Browne, 2001; Hendrickson & White, 1964), are available for factor rotation in the literature. In our simulation study, we applied promax rotation as it is one of the most computationally efficient oblique rotation methods in large-scale factor analysis. For the estimation implemented in the ‘mirt’ package, we use the built-in promax rotation to obtain the estimation, and for the pGVEM estimation, we use the function promax in R, with default $m=4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m=4$$\end{document} , to perform the rotation after the iteration ends.

The manipulated conditions include: (1) sample size $N=200,500$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=200,500$$\end{document} ; (2) test length $J=10,20$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=10,20$$\end{document} ; (3) number of categories $K=3,6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3,6$$\end{document} ; (5) low and high correlation among the latent traits; (6) small- and large-scale loadings. For each condition, a total number of 100 replicated cases were simulated. In the context of partial credit models, it is noted that the scaling of loadings plays a pivotal role in shaping the likelihood function. Specifically, when loadings are high and there exist multiple categories for partial credit scoring, the following case usually occurs: The probability of attaining the highest or lowest scores becomes disproportionately large. Consequently, this dominance of extreme scores may result in insufficient records of intermediate scores, thereby making the estimation of threshold parameters problematic. Therefore in the simulation studies, we considered two cases: (1) low scale loading: Parameter $a_{\mathrm{jr}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jr}$$\end{document} was simulated from Unif(0.5, 1) for all $j=1,\dots ,J,r=1,\dots ,D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,J,\;r=1,\ldots , D$$\end{document} ; (2) high scale loading: Parameter $a_{\mathrm{jr}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jr}$$\end{document} was simulated from Unif(1, 2) for all $j=1,\dots ,J,r=1,\dots ,D$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j=1,\ldots ,J,\;r=1,\ldots , D$$\end{document} . The threshold parameters $b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} are simulated from N(0, 1) for all $j=1,\dots ,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,\ldots , J$$\end{document} and $k=1,\dots ,K-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=1,\ldots , K-1$$\end{document} . For the latent variables, they were simulated from a multivariate normal distribution with 0 mean and 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 }}_\theta $$\end{document} . The diagonal entries were fixed as 1, and off-diagonal entries were generated from a uniform distribution Unif(0.1, 0.3) in the low correlation case and Unif(0.5, 0.7) in the high correlation case. For the responses generated from the simulated model parameters, we perform simulation only on cases where for all $j=1,\dots ,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,\ldots ,J$$\end{document} and $k=0,\cdots K-1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0,\cdots K-1$$\end{document} ,

$\begin{array}{c}\#\{i\mid Y_{\mathrm{ij}}=k,i=1,\dots ,N\}>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} \#\{i \mid Y_{ij}=k,i=1,\ldots ,N\}>0. \end{aligned}$$\end{document}

Here # denotes the set counting operator. Skipping any item with all responses being 0 is necessary since the threshold parameter linked to each category of this item cannot be identified within the finite sample context. For the convergence criterion, the algorithm was terminated when the change of all item parameters between two iterations dropped below a pre-specified threshold, i.e.,

(3.1) $\begin{array}{c}\frac{1}{J\times D+J\times K}\sum _{j=1}^J\left(\sum _{r=1}^D,{\left(a_{\mathrm{jr}}^{\left(t\right)},-,a_{\mathrm{jr}}^{(t-1)}\right)}^2,+,\sum _{k=1}^{K-1},{\left(b_{\mathrm{jk}}^{\left(t\right)},-,b_{\mathrm{jk}}^{(t-1)}\right)}^2\right)<{10}^{-5}.\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}{J\times D+J\times K}\sum _{j=1}^J\left[ \sum _{r=1}^D\left( a_{jr}^{(t)}- a_{jr}^{(t-1)}\right) ^2+\sum _{k=1}^{K-1}\left( b_{jk}^{(t)}-b_{jk}^{(t-1)}\right) ^2\right] <10^{-5}. \end{aligned}$$\end{document}

The estimation errors are presented separately in the form of mean squared error and bias for the discrimination and threshold parameters and the covariance matrix, averaged across the test items:

$\begin{array}{cc}Bia{s}_a=\frac{1}{\mathrm{JD}}\sum _{j=1}^J\sum _{r=1}^D{\hat{a}}_{\mathrm{jr}}-a_{\mathrm{jr}},& MS{E}_a=\frac{1}{\mathrm{JD}}\sum _{j=1}^J{\Vert a_j-{\hat{a}}_j\Vert }_2^2;\\ Bia{s}_b=\frac{1}{J(K-1)}\sum _{j=1}^J\sum _{k=1}^{K-1}{\hat{b}}_{\mathrm{jr}}-b_{\mathrm{jr}},& MS{E}_b=\frac{1}{J(K-1)}\sum _{j=1}^J\sum _{k=1}^{K-1}{(b_{\mathrm{jk}}-{\hat{b}}_{\mathrm{jk}})}^2;\\ Bia{s}_{\Sigma }=\frac{2}{D(D-1)}\sum _{l<h}{\hat{\Sigma }}_{\mathrm{hl}}-{\Sigma }_{\mathrm{hl}},& MS{E}_{\Sigma }=\frac{2}{D(D-1)}{\Vert {\Sigma }_{\theta }-{\hat{\Sigma }}_{\theta }\Vert }_F^2.\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} Bias_a=\frac{1}{JD}\sum _{j=1}^{J}\sum _{r=1}^{D} \hat{ a}_{jr}-a_{jr},\quad&MSE_a=\frac{1}{JD}\sum _{j=1}^{J}\Vert {\varvec{a}}_j-\hat{{\varvec{a}}}_j\Vert ^2_2;\\ Bias_b =\frac{1}{J(K-1)}\sum _{j=1}^{J}\sum _{k=1}^{K-1} \hat{ b}_{jr}-b_{jr},\quad&MSE_b=\frac{1}{J(K-1)}\sum _{j=1}^J\sum _{k=1}^{K-1}(b_{jk}-{\hat{b}}_{jk})^2;\\ Bias_\Sigma =\frac{2}{D(D-1)}\sum _{l<h}{\hat{\Sigma }}_{hl}-\Sigma _{hl},\quad&MSE_\Sigma =\frac{2}{D(D-1)}\Vert {\varvec{\Sigma }}_\theta -\hat{{\varvec{\Sigma }}}_\theta \Vert _F^2. \end{aligned}$$\end{document}

Note that here both the true and estimated discrimination parameters have been rotated by the promax rotations. The number of Markov chain samples drawn in the MHRM algorithm was by default 5, 000 in the R package “mirt.” The convergence criterion and optimizer were all set as default in the ‘mirt’ package.

Figure 1 Mean squared error of the estimation for the multidimensional generalized partial credit model of 3 categories from exploratory factor analysis with small-scale loadings using different methods.

Figure 2 Mean squared error of the estimation for multidimensional generalized partial credit model of 6 categories from exploratory factor analysis with small-scale loadings using different methods.

Figure 3 Bias of the estimation for the multidimensional generalized partial credit model of 3 categories from exploratory factor analysis with small-scale loadings using different methods.

Figure 4 Bias of the estimation for the multidimensional generalized partial credit model of 6 categories from exploratory factor analysis with small-scale loadings using different methods.

Under the setting of small-scale loadings (i.e., $a_{\mathrm{jr}}\sim$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jr}\sim $$\end{document} Unif(0.5, 1)), Figs. 1 and 2 show the MSE of the four estimation methods for 3-category and 6-category cases, respectively. Each box represents the distribution of errors from 100 replications. We truncated the scale of the y-axis of the plot to make it easier to compare the estimation precision across different scenarios. We present the full version of the boxplots in Appendix B. Overall, our method provides more stable and accurate estimates of the MGPCM, as seen from both figures. The observed results indicate a reduced variability in estimation errors for the model parameters when comparing pGVEM to the other methods. Both MHRM and StEM exhibit better stability and accuracy compared to the standard EM algorithm. Notably, pGVEM demonstrates good stability, particularly evident when the sample size N is 200. Additionally, it is noteworthy that as the number of categories increases, the estimation of the threshold parameters becomes more challenging, an anticipated result given that the model becomes more complicated for multicategory cases. Interestingly, despite a decrease in accuracy, the proposed pGVEM method exhibits a comparatively modest increase in variability in many cases compared with alternative methods, indicating the capability of the pGVEM method to handle more complex scenarios. Furthermore, we present the bias of the estimation using the four different methods in Figs. 3 and 4 for the considered cases. In general, the bias observed in pGVEM estimation tends to be more moderate across various cases, particularly with regard to the threshold parameter.

Figure 5 Mean squared error of the estimation for the multidimensional generalized partial credit model of 3 categories from exploratory factor analysis with large-scale loadings using different methods.

Figure 6 Mean squared error of the estimation for multidimensional generalized partial credit model of 6 categories from exploratory factor analysis with large-scale loadings using different methods.

Figure 7 Bias of the estimation for the multidimensional generalized partial credit model of 3 categories from exploratory factor analysis with large-scale loadings using different methods.

Figure 8 Bias of the estimation for the multidimensional generalized partial credit model of 6 categories from exploratory factor analysis with large-scale loadings using different methods.

Figure 9 Computational time comparison for four methods with $K=6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=6$$\end{document} : left side depicts $J=10$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=10$$\end{document} , while the right side corresponds to $J=20$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=20$$\end{document} .

Under the setting of large-scale loadings (i.e., $a_{\mathrm{jr}}\sim$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jr}\sim $$\end{document} Unif(1,2)), the MSE and bias results are presented in Figs. 5, 6, 7, and 8. We can see from the simulation results that in this setting the error of estimation gets larger in most cases. The reason is that high level of loading increases the frequency of extreme scores, thus making the estimation more challenging. Yet pGVEM still outperforms the other methods with the regime of small sample size or low correlation. Furthermore, in terms of the recovery of the threshold parameter, pGVEM provides estimates with less bias and error. It is noteworthy that in this scenario, it occurs more frequently that for some item, there is no record of certain category from any individual, leading to us discarding such cases from the analysis. This suggests that when dealing with multiple categories of responses, the discrimination parameters may tend to be small for model interpretability.

In Fig. 9 we present the averaged computation time of the four methods under different settings of sample size N and test length J. The results exhibit similarity across the cases, and for brevity, we displayed the case where the discrimination parameters are simulated from Unif(0.5, 1) and correlation coefficient from Unif(0.1, 0.3). We take a total number of 6 categories. It is obvious that, compared with pGVEM, the traditional EM algorithm is inefficient especially when the sample size is large. Also, pGVEM is slightly faster than the stochastic EM algorithm and achieves similar computation efficiency compared with MH-RM algorithm in the displayed case. The computational efficiency of our pGVEM algorithm makes it possible to provide fast estimation on large datasets.

The results of our study demonstrate the superiority of the pGVEM algorithm over the traditional EM algorithm along with MH-RM and StEM in terms of parameter recovery and computational efficiency. Specifically, the pGVEM algorithm achieves comparable, and in some cases, superior parameter recovery compared to the other algorithms. Moreover, pGVEM generates fewer estimation outliers, particularly in situations where the sample size and test length are large.

As a side check, we compared the pGVEM algorithm with the GVEM algorithm proposed by Cho et al. (2021) for the special case of M2PL. The comparison was made under identical setting in this section, except for that we take $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 results, presented in Appendix B, indicate that our algorithm performs similarly to Cho et al. (2021) with binary data. Overall, our findings demonstrate that the pGVEM algorithm is a robust and efficient method for parameter recovery for the MGPCM. It can also be applied to other models, including M2PL model, with similar success. The results suggest that the pGVEM algorithm may be a valuable tool for researchers seeking to analyze complex data structures efficiently and accurately.

3.2. Study II

We conducted simulation study II to assess the performance of the proposed bootstrap standard error (SE) estimation procedure. We explore the bootstrap estimation in the multidimensional case with multiple categories. In the multidimensional case, as clarified in Sect. 2.2, we found that the traditional methods (EM and MH-RM) exhibited instability and produced infeasible results across numerous settings. Consequently, in this section, we focus on the bootstrap-based SE estimates under the EM algorithm, MH-RM, and the proposed pGVEM algorithm.

The comparisons were conducted under the simulation setting similar to Simulation Study I and the manipulated factors include sample size, test length, factor correlations, and the number of categories. The empirical standard deviations of the estimated item parameters as

$\begin{array}{c}S{E}_v=\frac{1}{R-1}\sum _{r=1}^R{({\hat{v}}^{\left(r\right)}-v)}^2\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} SE_{v}=\frac{1}{R-1}\sum _{r=1}^R({\hat{v}}^{(r)}-v)^2 \end{aligned}$$\end{document}

across replications per condition are considered as the approximations of true SEs for each method. Here v stands for item parameters to represent $a_{\mathrm{jd}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jd}$$\end{document} or $b_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} , and ${\hat{v}}^{\left(r\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{v}}^{(r)}$$\end{document} is the estimated parameter in the rth replicate. To assess the performance of the proposed method, we computed SE estimations along with their bias and relative bias as follows:

$\begin{array}{cc}AverageSE& =\frac{1}{J(D+K)}\sum _{j=1}^J\left(\sum _{r=1}^D,{\widehat{\mathrm{SE}}}_{a_{\mathrm{jr}}},+,\sum _{k=1}^{K-1},{\widehat{\mathrm{SE}}}_{b_{\mathrm{jk}}}\right)\\ Bias& =\frac{1}{J(D+K)}\sum _{j=1}^J\left(\sum _{r=1}^D,{\widehat{\mathrm{SE}}}_{a_{\mathrm{jr}}},-,S,E_{a_{\mathrm{jr}}},+,\sum _{k=1}^{K-1},{\widehat{\mathrm{SE}}}_{b_{\mathrm{jk}}},-,S,E_{b_{\mathrm{jk}}}\right)\\ RelativeBias& =\frac{1}{J(D+K)}\sum _{j=1}^J\left(\sum _{r=1}^D,({\widehat{\mathrm{SE}}}_{a_{\mathrm{jr}}}-S{E}_{a_{\mathrm{jr}}}),/,S,E_{a_{\mathrm{jr}}},+,\sum _{k=1}^{K-1},({\widehat{\mathrm{SE}}}_{b_{\mathrm{jk}}}-S{E}_{b_{\mathrm{jk}}}),/,S,E_{b_{\mathrm{jk}}}\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} Average \;SE&=\frac{1}{J(D+K)}\sum _{j=1}^J\left[ \sum _{r=1}^D{\widehat{SE}}_{a_{jr}}+\sum _{k=1}^{K-1}{\widehat{SE}}_{b_{jk}}\right] \\ Bias&=\frac{1}{J(D+K)}\sum _{j=1}^J\left[ \sum _{r=1}^D{\widehat{SE}}_{a_{jr}}-SE_{a_{jr}}+\sum _{k=1}^{K-1}{\widehat{SE}}_{b_{jk}}-SE_{b_{jk}}\right] \\ Relative\;Bias&=\frac{1}{J(D+K)}\sum _{j=1}^J\left[ \sum _{r=1}^D({\widehat{SE}}_{a_{jr}}{-}SE_{a_{jr}})/SE_{a_{jr}}{+}\sum _{k=1}^{K-1}({\widehat{SE}}_{b_{jk}}{-}SE_{b_{jk}})/SE_{b_{jk}}\right] , \end{aligned}$$\end{document}

providing a comprehensive assessment of the reliability and accuracy of the proposed bootstrap methods. Here we present the SEs of $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} pooled together with the aim of showing the effectiveness of our method compared with its alternatives. The results 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} are similar to the pooled ones. The variability of discrimination and threshold parameters has been well illustrated in the previous study.

The results are shown in Fig. 10 for the low correlation setting and 11 for the high correlation setting. In the multidimensional case, where most of the methods fail to estimate the SEs numerically or from Fisher’s information matrix, the bootstrap method still generates stable results. In comparison with bootstrap methods from alternative estimations, the pGVEM-based bootstrap exhibits a lower bias. It is observed that when the sample size is 200, the pGVEM-based bootstrap may slightly underestimate standard errors. Nevertheless, its overall performance remains relatively strong across diverse settings.

Figure 10 Standard error assessment of estimation for the multidimensional generalized partial credit model from exploratory factor analysis using different methods with $D=3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=3$$\end{document} in low factor correlation setting.

Figure 11 SE assessment of estimation for the multidimensional generalized partial credit model from exploratory factor analysis using different methods with $D=3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=3$$\end{document} in high factor correlation setting.

3.3. Study III

In this section we conduct a simulation study to examine the performance of the proposed AIC and BIC in determining the latent dimension. Considering the complexity of the model setting of MGPCM, we explore the accuracy of factor identification in sample size $N=500,800$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=500,800$$\end{document} and test length $J=20,40$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=20,40$$\end{document} . We consider both low correlation settings (simulated from Unif(0.1, 0.3)) and high correlation settings (simulated from Unif(0.5, 0.7)) for dimensions $D=2,3,4$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=2,3,4$$\end{document} . In each configuration, a total number of 100 independent samples were generated and we recorded the number of cases where the number of factors was correctly identified. Discrimination parameters are generated from Unif(1, 2).

Tables 1 and 2 present the correct estimation rates for the number of dimensions of 3 and 6 categories, respectively. The results indicate that an increase in sample size generally contributes to higher correct estimation rates. We can also observe that a lower correlation generally leads to higher correct estimation rates, especially when there are fewer categories and shorter test lengths, as the latent structure may be more challenging to identify in such settings. Overall, we observe that AIC performs slightly better in the case of $K=3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document} and BIC has an advantage in the regime of $K=6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=6$$\end{document} . In conclusion, the proposed criterion is efficient in general, and with a larger sample size, the model is more likely to be correctly identified. Our findings are also consistent with existing studies under the MIRT model (Cho et al., 2021).

4. Real Data Analysis

4.1. Trend in International Mathematics and Science Study Dataset

In this section, we demonstrate the application of the pGVEM algorithm by analyzing a dataset from the Trend in International Mathematics and Science Study (TIMSS) (Mullis & Martin, 2017; Martin & Mullis, 2019; Fishbein et al., 2018; Martin et al., 2020). TIMSS provides reliable and timely trend data on the mathematics and science achievement of US students compared to that of students in other countries. The assessment consists of a large pool of mathematics and science questions, which are divided into different blocks using the item matrix sampling design to relieve response burden. Specifically, 14 booklets were assembled, and each student was required to complete one of them. Different booklets may include the same items for linking. TIMSS 2019 divided the test items into 28 blocks for each grade, with 13 mathematics items and 15 science items. In the data of students in grade eight, each block contained 12 to 18 items. For our analysis, we selected a mathematics and a science block that appeared in booklet 5 of grade 8. In the mathematics block there were 2 polytomous items, and in the science there were 3. Of the 1,252 students who responded to these booklets, 918 students’ responses were completely recorded. In this study, we only estimated the parameters using the data from these 918 students. Appendix C provides details of the item code and corresponding test content. The IRT parameters provided in the TIMSS assessment document (Martin & Mullis, 2019) were used as the true parameters, to which our estimated parameters were compared.

It should be noted that in operational analysis of TIMSS, uni-dimensional IRT models were used for math and science domains separately. When analyzing the items from math and science domain separately, the modified information criterion in Equation (2.14) and (2.15) computed under the EFA framework both attain the smallest value when the latent dimension is 1 (i.e., $D=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=1$$\end{document} ), which implies that both domains are essentially uni-dimensional under the generalized partial credit model setting. However, according to the information criterion provided in the ‘mirt’ package (i.e., using its default EM algorithm), the latent dimension cannot be clearly decided. In the following we display the parameters estimated by EM with fixed quadrature and pGVEM. The results are as follows.

To show the results visually, we plot the estimated parameters, by EM and pGVEM, in Figs. 12 and 13. Based on our analysis, it can be inferred that the b-parameters obtained from both methods exhibit a significant level of proximity. Similarly, the a-parameters demonstrate close values. However, it is worth noting that the estimates derived from the EM algorithm tend to be slightly larger for the majority of the items.

Table 1 Correct number of trials in determining the latent dimension, $K=3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document}

Table 2 Correct number of trials in determining the latent dimension, $K=6$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=6$$\end{document}

Figure 12 Comparison of discrimination parameters estimated from ‘mirt’ and pGVEM.

To reveal the relation of latent ability across different domains, we estimate all 28 items jointly. First, compute the modified information criterion in Equation (2.14) and (2.15) for the setting of joint estimating the data from mathematics and science block, with results presented in Table 3 and 4. From the results the optimal number of latent dimensions is 2. This result is reasonable as math proficiency and science proficiency should emerge as two separate factors. Yet the information criterion output from ‘mirt’ package again did not provide sensible results. By assuming the latent dimension is 2, we present the results estimated by the two methods after promax rotation in Fig. 14. The analysis of estimated parameters reveals a clear two-factor loading structure as shown in Fig. 14. Notably, the second dimension emerges as primarily the students’ proficiency in solving mathematical problems, whereas the first dimension tends to measure students’ science proficiency. This finding suggests that the two dimensions capture distinct but interconnected aspects of overall cognitive ability. Moreover, we can infer that certain types of questions are particularly effective in assessing students’ latent abilities in either mathematics or science. These questions demonstrate a higher degree of typicality in evaluating students’ competence within their respective domains. This insight underscores the importance of carefully selecting assessment items that align with the targeted cognitive skills, as they offer a more accurate reflection of students’ underlying abilities.

Figure 13 Comparison of threshold parameters estimated from ‘mirt’ and pGVEM.

Table 3 Information criterion from ‘mirt’. The smallest value of each information criterion is given in bold

Math, Sci, Joint denote, respectively, information criterion for math items, science items and the collection of all items jointly.

Table 4 Information criterion from pGVEM algorithm. The smallest value of each information criterion is given in bold

Math, Sci, Joint denote, respectively, information criterion for math items, science items and the collection of all items jointly.

Figure 14 Jointly estimated discrimination parameters from ‘mirt’ and pGVEM.

4.2. Big-Five Factor Personality Dataset

We further demonstrate the application of the pGVEM algorithm by analyzing a dataset from the Big-Five personality assessment. The five-factor model (FFM) stands as the most widely recognized and extensively used model in psychology for understanding and measuring personality (Goldberg, 1992). This model provides a framework to capture the richness and complexity of individual differences, assuming five domains that encompass a comprehensive spectrum of personality traits. The five dimensions include Openness, Conscientiousness, Extraversion, Agreeableness, and Neuroticism (often referred to by the acronym OCEAN) (Costa & McCrae, 2008).

To measure the latent traits, various assessment tools have been developed to assess an individual’s standing on the Big-Five dimensions (Wiggins & Trapnell, 1997; McCrae et al., 1996; Goldberg et al., 1999). In our study, we employ the Big-Five Factor Markers derived from the International Personality Item Pool (IPIP), a widely recognized instrument developed by Goldberg (1992). The dataset is publicly available at http://openpsychometrics.org/tests/IPIP-BFFM/. This dataset consists of responses from a substantial sample of 19,718 individuals, each evaluated on the fifty-item Big-Five Factor Markers. The IPIP consists of fifty items, each requiring respondents to rate the extent to which they perceive each statement as true about themselves on a five-point scale. This scale ranges from 1 (Disagree) to 3 (Neutral) to 5 (Agree), providing a nuanced and graded assessment across the five categories.

Under the setting of exploratory factor analysis, we fit the MGPCM using the proposed pGVEM algorithm. The estimated factor loadings are represented in Fig. 15. The heatmap reveals that each factor exhibits a distinct and salient association with the grouped items. Also, we can see that not all latent factors have a positive influence on the overall rating—a phenomenon commonly observed in personality tests (McCrae et al., 1996; Goldberg et al., 1999). A significant correlation between Factor 1 and Factor 3 can be observed, as they both significantly impact responses to the A-items. To further illustrate these relationships, we present the estimated correlation matrix of the factors in Table 5. The largest correlation is between Factor 1 and Factor 3 with a value of 0.529, which is consistent with our observation from Fig. 15. Overall, the estimated factor structure aligns with existing literature, demonstrating the usefulness of the proposed method as a computationally efficient tool to analyze large scale assessment data.

Figure 15 Heatmap of loading matrix from Big-Five personality study.

Table 5 Estimated correlation matrix of the 5 latent factors for the Big-Five dataset. The absolute values of correlations higher than 0.1 are given in bold

5. Conclusion

In this paper, we proposed a new Gausssian variational estimation algorithm of polytomous responses (namely, pGVEM) for the multidimensional generalized partial credit model (MGPCM). The MGPCM is one of the most widely used models in cognitive assessment and educational measurement for items that are scored polytomously, such as assigning partial scores for intermediate correct steps. There are limited methods for an efficient estimation of MGPCM, and it is shown that our pGVEM algorithm outperforms many existing approaches in that it is relatively stable and much faster without resulting in much aberrant estimates or convergence issues. When the sample size and test length are both large, our proposed variational lower bound seems to approximate the target marginal likelihood closely. The computation efficiency is achieved by replacing the intractable high-dimensional integral with a variational lower bound that contributes to faster EM-type updates involving only small-scale linear equations. The simulation study in Sect. 3 provides simulation evidence to support pGVEM in producing accurate parameter estimates quickly, as compared to the traditional EM implementation of marginal maximum likelihood estimators. The real data analysis demonstrates that our estimation scheme is capable of extracting proper information about the latent variables.

Variational inference has emerged as a prominent and efficient methodology in the field of psychometrics, particularly due to its ability to handle large-scale datasets with both accuracy and computational efficiency. In addition to IRT models with continuous latent variables, Yamaguchi and Okada (2020) recently proposed a variational Bayesian (VB) inference algorithm for the saturated cognitive diagnosis models, which represents a notable advancement in scalable and computationally efficient Bayesian estimation for discrete latent variable models. Oka and Okada (2023) developed a scalable estimation algorithm for the DINA Q-matrix, which employs an iteration scheme utilizing stochastic optimization and variational inference. Our method, on the other hand, extends the literature on continuous latent variable estimation (Cho et al., 2021, 2022) by considering multiple response categories. It is encouraging to further explore the relationship of various variational approximation methods, which may lead to a more robust and flexible estimation framework for many psychometric models.

There are also some other potential directions to extend the current work. First, as with many nonconvex optimization problems, the efficacy of our algorithm may be influenced by the chosen initial values. Instances where initial values deviate substantially from the actual parameters may lead the algorithm to converge to local optima rather than the global one. Given the propensity for variational approximations to yield multiple local optimal values, investigating the initialization strategy for our estimation process warrants further exploration. Second, the current method of determining the number of latent dimensions may not work for certain cases, especially when the dimension is very high or there exists a high correlation between latent factors. Therefore, a more accurate and robust model selection method that will work in these challenging scenarios is needed.