Additional Penalty for M3PL

The parameter estimation for M3PL in practice often gets more challenging due to the inclusion of guessing parameters. To tackle this challenge and improve the accuracy of the parameter estimation in M3PL, we propose to impose additional constraints on the model parameters, $B=\{b_j;j=1,\cdots ,J\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {B}}= \{b_j; j=1, \dots , J\}$$\end{document} and $C=\{c_j;j=1,\cdots ,J\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {C}}= \{c_j; j=1,\dots , J\}$$\end{document} in addition to the parameter matrix $A$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {A}}$$\end{document} . Specifically for parameter estimation, we solve the following optimization problem where $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 )$$\end{document} denotes a penalty function on model parameters:

(6) $\begin{array}{c}({\hat{A}}_{\lambda },{\hat{B}}_{\lambda },{\hat{C}}_{\lambda })={\text{argmax}}_{A,B,C}l(A,B,C;Y)-P_{\lambda }\left(A\right)+P\left(B\right)+P\left(C\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} ({\hat{{\mathbf {A}}}}_{\lambda }, {\hat{{\mathbf {B}}}}_{\lambda }, {\hat{{\mathbf {C}}}}_{\lambda }) = \text {argmax}_{{\mathbf {A}},{\mathbf {B}},{\mathbf {C}}} \qquad {l({\mathbf {A}},{\mathbf {B}},{\mathbf {C}}; {\mathbf {Y}}) - P_\lambda ({\mathbf {A}}) + P({\mathbf {B}}) + P({\mathbf {C}})} \end{aligned}$$\end{document}

where $P\left(B\right)={\sum }_{j=1}^{J}logN\left(b_j\right|{\mu }_b,{\sigma }_b^2)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\mathbf {B}}) = \sum _{j=1}^J\log N(b_j|\mu _b,\sigma _b^2)$$\end{document} , and $P\left(C\right)={\sum }_{j=1}^{J}logBeta\left(c_j\right|{\alpha }_c,{\beta }_c)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P({\mathbf {C}}) = \sum _{j=1}^J\log Beta(c_j|\alpha _c, \beta _c)$$\end{document} for some distribution parameters ${\mu }_b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _b$$\end{document} , ${\sigma }_b^2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _b^2$$\end{document} , ${\alpha }_c$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _c$$\end{document} , and ${\beta }_c$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _c$$\end{document} . These penalty functions are chosen to satisfy the ranges of values on which the parameters are defined. For instance, since the guessing parameters $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {C}}$$\end{document} naturally satisfy the constraint $\{0<c_j<1;j=1,\cdots ,J\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{0< c_j < 1; j=1,\dots ,J\}$$\end{document} , we can assume a “prior” distribution of $c_j\sim Beta({\alpha }_c,{\beta }_c)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_j \sim Beta(\alpha _c,\beta _c)$$\end{document} . Similarly, we can assume a “prior” distribution of $b_j\sim N({\mu }_b,{\sigma }_b^2)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j \sim N(\mu _b,\sigma ^2_b)$$\end{document} . The penalty on $b_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} and $c_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_j$$\end{document} are essentially a $L_2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} -type and Laplace penalization, respectively. By imposing these additional penalties on model parameters $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {B}}$$\end{document} and $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {C}}$$\end{document} , the parameter estimation becomes more stable and robust.

The approach of imposing additional penalty on model parameters $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {B}}$$\end{document} and $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {C}}$$\end{document} with the chosen distributions is similar to the Bayes modal estimation presented by Tierney and Kadane (Reference Tierney and Kadane1986). That is, an augmented optimization objective is employed that includes the likelihood and some prior beliefs on the item parameters. These priors can be used to prevent deviant parameter estimates and help the algorithm to produce more accurate estimation in complex M3PL models. Essentially, Bayes modal estimation can be seen as a regularization on maximum likelihood estimation where maximum likelihood estimation is a special case of Bayes model estimation that assumes uniform prior distributions.

The amount of penalization can be flexibly controlled using the distribution parameters. For instance, one can use non-informative priors on $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {C}}$$\end{document} such as Beta(1, 1), which is equivalent to flat uniform distribution on [0, 1]. Additionally, one can similarly choose non-informative normal prior with high variance ${\sigma }_b$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _b$$\end{document} for $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {B}}$$\end{document} . This suggests that although additional penalization functions are added, the algorithm also allows the flexible estimation with essentially no penalty with the choice of non-informative distributions. The advantage of this is that practitioners can adjust the amount of prior knowledge they would like to impose on the model. The less prior knowledge one uses, the more flexible the estimation is and the results will be based more on the observed data. With these prior-like penalties, our algorithm yields more precise parameter estimates for the M3PL model.

Computation via GVEM

This section introduces the main estimation algorithm to obtain the estimate $({\hat{A}}_{\lambda },{\hat{B}}_{\lambda },{\hat{C}}_{\lambda })$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\hat{A}}_{\lambda }, {\hat{B}}_{\lambda }, {\hat{C}}_{\lambda })$$\end{document} via (6) using GVEM algorithm. As introduced in Sect. 2, we will use a variational lower bound to approximate the intractable marginal log-likelihood $l(A,B,C;Y)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l({\mathbf {A}},{\mathbf {B}},{\mathbf {C}}; {\mathbf {Y}})$$\end{document} in (6).

To derive a lower bound for easy estimation of the M3PL parameters, instead of directly working with (4), we employ an equivalent representation of the M3PL model with auxiliary latent variable $Z_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij}$$\end{document} , which is an indicator function of whether the ith individual answers the jth item based on the latent ability or guesses it correctly (von Davier, Reference von Davier2009). Specifically $Z_{\mathrm{ij}}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij}=1$$\end{document} if the ith individual solves item j based on his/her ability, and $Z_{\mathrm{ij}}=0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij} = 0$$\end{document} if he/she guesses item j correctly. The distribution of $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} given 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} and $Z_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij}$$\end{document} is then

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

where we define $0^0=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0^0=1$$\end{document} , and it can be seen that this new model with auxiliary variable Z is equivalent to the M3PL model (von Davier, Reference von Davier2009, Cho et al., Reference Cho, Wang, Zhang and Xu2021). Denote $Z_i=\{Z_{i1},Z_{i2},\cdots ,Z_{\mathrm{iJ}}\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{Z}}_i = \{Z_{i1}, Z_{i2}, \dots , Z_{iJ} \}$$\end{document} and its distribution as $p\left(Z_i\right)={\prod }_{j=1}^{J}p\left(Z_{\mathrm{ij}}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p({\varvec{Z}}_i) = \prod _{j=1}^J p(Z_{ij})$$\end{document} . Then the complete data likelihood of the ith subject can be written as

(7) $\begin{array}{cc}& logP(Y_i,{\theta }_i,Z_i\mid A,B,C)\\ =& logP(Y_i\mid {\theta }_i,Z_i,A,B,C)+log\varphi \left({\theta }_i\right)+logp\left(Z_i\right)\\ =& \sum _{j=1}^J\left(Y_{\mathrm{ij}},Z_{\mathrm{ij}},({\alpha }_j^{\top }{\theta }_i-b_j),+,Z_{\mathrm{ij}},log,\frac{1}{1+exp({\alpha }_j^{\top }{\theta }_i-b_j)},+,(1-Z_{\mathrm{ij}}),log,I,(Y_{\mathrm{ij}}=1)\right)\\ & +log\varphi \left({\theta }_i\right)+logp\left(Z_i\right),\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\log P(Y_i,{\varvec{\theta }}_i, {\varvec{Z}}_i \mid {\varvec{A}}, {\varvec{B}}, {\varvec{C}})\nonumber \\= & {} \log P(Y_i\mid {\varvec{\theta }}_i, {\varvec{Z}}_i, {\varvec{A}},{\varvec{B}},{\varvec{C}}) +\log \phi ({\varvec{\theta }}_i) + \log p({\varvec{Z}}_i)\nonumber \\= & {} \sum _{j=1}^J \left\{ Y_{ij} Z_{ij} ({\varvec{\alpha }}_j^\top {\varvec{\theta }}_i-b_j) + Z_{ij} \log \frac{1}{1+\exp ({\varvec{\alpha }}_j^\top {\varvec{\theta }}_i-b_j)} + (1-Z_{ij})\log I(Y_{ij}=1)\right\} \nonumber \\&+ \log \phi ({\varvec{\theta }}_i) + \log p({\varvec{Z}}_i), \end{aligned}$$\end{document}

where $\varphi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} denotes the normal probability density function for latent variable $\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} . Here, without loss of generality, we focus on the ith subject’s likelihood function due to the independence of different subjects.

With the above representation, for any variational distribution functions $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$$\end{document} and $r_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}$$\end{document} (to be estimated later) of 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} and $Z_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_{ij}$$\end{document} , similar to the derivation in Sect. 2, we have the following variational lower bound, which generalizes (4),

(8) $\begin{array}{cc}logP(Y_i\mid A,B,C)\ge & {\int }_{{\theta }_i}\sum _{Z_i}logP(Y_i,{\theta }_i,Z_i\mid A,B,C)\times q_i\left({\theta }_i\right)r_i\left(Z_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} \log P(Y_{i}\mid {\varvec{A}},{\varvec{B}},{\varvec{C}})\ge & {} \int _{{\varvec{\theta }}_i} \sum _{{\varvec{Z}}_i}\log P(Y_{i},{\varvec{\theta }}_i, {\varvec{Z}}_i\mid {\varvec{A}},{\varvec{B}},{\varvec{C}}) \times q_i({\varvec{\theta }}_i)r_i({\varvec{Z}}_i)d{\varvec{\theta }}_i \end{aligned}$$\end{document}

(9) $\begin{array}{cc}& -{\int }_{{\theta }_i}\sum _{Z_i}log(q_i\left({\theta }_i\right)r_i\left(Z_i\right))\times q_i\left({\theta }_i\right)r_i\left(Z_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}&- \int _{{\varvec{\theta }}_i} \sum _{{\varvec{Z}}_i} \log \big (q_i({\varvec{\theta }}_i)r_i({\varvec{Z}}_i)\big )\times q_i({\varvec{\theta }}_i)r_i({\varvec{Z}}_i)d{\varvec{\theta }}_i, \end{aligned}$$\end{document}

where $r_i\left(Z_i\right)={\prod }_{j=1}^J{r}_{\mathrm{ij}}\left(Z_{\mathrm{ij}}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_i({\varvec{Z}}_i) = \prod _{j=1}^J r_{ij}(Z_{ij})$$\end{document} . Since (9) doesn’t depend on parameters $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} , $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{B}}$$\end{document} and $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{C}}$$\end{document} , we focus on (8) for the derivation of the lower bound. For (8), note that $logP(Y_i,{\theta }_i,Z_i\mid A,B,C)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log P(Y_{i},{\varvec{\theta }}_i, {\varvec{Z}}_i\mid {\varvec{A}},{\varvec{B}},{\varvec{C}})$$\end{document} takes the form of (7). To obtain a closed form lower bound expression for (8), we further use a local variational method (Bishop, Reference Bishop2006; Jordan, Ghahramani, Jaakkola, & Saul, Reference Jordan, Ghahramani, Jaakkola and Saul1999). Particularly, define ${\xi }_{i,j}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _{i,j}$$\end{document} as a variational parameter indexed by i and j, and let $\eta \left({\xi }_{i,j}\right)={\left(2{\xi }_{i,j}\right)}^{-1}[e^{{\xi }_{i,j}}/(1+e^{{\xi }_{i,j}})-1/2].$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \eta ({\xi _{i,j}})= (2{\xi _{i,j}})^{-1}[ {e^{\xi _{i,j}}}/{(1+e^{\xi _{i,j}})}-{1}/{2} ]. $$\end{document} Let ${\xi }_i=({\xi }_{i,j},j=1,\cdots ,J)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\xi }_i = (\xi _{i,j}, j=1,\cdots , J)$$\end{document} denote the ith subject’s variational parameters for the J items. Then following the local variational method (Bishop, Reference Bishop2006), we have

$\begin{array}{cc}& logP(Y_i,{\theta }_i,Z_i\mid A,B,C)\ge l(A,B,C,{\xi }_i;Y_i,{\theta }_i,Z_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}&\log P(Y_i,{\varvec{\theta }}_i, {\varvec{Z}}_{i} \mid {\varvec{A}}, {\varvec{B}}, {\varvec{C}}) \ge l({\varvec{A}}, {\varvec{B}},{\varvec{C}}, \varvec{\xi }_i ; Y_i,{\varvec{\theta }}_i,{\varvec{Z}}_{i}), \end{aligned}$$\end{document}

where $l(A,B,C,{\xi }_i;Y_i,{\theta }_i,Z_i)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l({\varvec{A}}, {\varvec{B}},{\varvec{C}}, \varvec{\xi }_i ; Y_i,{\varvec{\theta }}_i,{\varvec{Z}}_{i})$$\end{document} is defined as

(10) $\begin{array}{cc}& l(A,B,C,{\xi }_i;Y_i,{\theta }_i,Z_i)\\ =& \sum _{j=1}^J{Z}_{\mathrm{ij}}log\frac{e^{{\xi }_{i,j}}}{(1+e^{{\xi }_{i,j}})}+\sum _{j=1}^J{Z}_{\mathrm{ij}}{Y}_{\mathrm{ij}}({\alpha }_j^{\top }{\theta }_i-b_j)+\sum _{j=1}^J\frac{1}{2}{Z}_{\mathrm{ij}}(b_j-{\alpha }_j^{\top }{\theta }_i-{\xi }_{i,j})\\ & -\sum _{j=1}^J{Z}_{\mathrm{ij}}\eta \left({\xi }_{i,j}\right)\{{(b_j-{\alpha }_j^{\top }{\theta }_i)}^2-{\xi }_{i,j}^2\}\\ & +\sum _{j=1}^J\left((1-Z_{\mathrm{ij}}),log,I,(Y_{\mathrm{ij}}=1)\right)+log\varphi \left({\theta }_i\right)+logp\left(Z_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}&l({\varvec{A}}, {\varvec{B}},{\varvec{C}}, \varvec{\xi }_i ; Y_i,{\varvec{\theta }}_i,{\varvec{Z}}_{i}) \nonumber \\= & {} \sum _{j=1}^J Z_{ij}\log \frac{e^{\xi _{i,j}}}{(1+e^{\xi _{i,j}})} +\sum _{j=1}^J Z_{ij} Y_{ij}({\varvec{\alpha }}_j^\top {\varvec{\theta }}_i-b_j) +\sum _{j=1}^J \frac{1}{2}Z_{ij}(b_j-{\varvec{\alpha }}_j^\top {\varvec{\theta }}_i-{\xi _{i,j}})\nonumber \\&-\sum _{j=1}^J Z_{ij}\eta ({\xi _{i,j}})\{(b_j-{\varvec{\alpha }}_j^\top {\varvec{\theta }}_i)^2-\xi _{i,j}^2\}\nonumber \\&+ \sum _{j=1}^J \left\{ (1-Z_{ij})\log I(Y_{ij}=1) \right\} +\log \phi ({\varvec{\theta }}_i) + \log p({\varvec{Z}}_{i}), \end{aligned}$$\end{document}

and it gives a lower bound of $logP(Y_i,{\theta }_i,Z_i\mid A,B,C)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log P(Y_i,{\varvec{\theta }}_i, {\varvec{Z}}_{i} \mid {\varvec{A}}, {\varvec{B}}, {\varvec{C}})$$\end{document} in (8). We then have the following expression for the variational lower bound of the marginal likelihood of all observed responses in (6),

$\begin{array}{cc}l(A,B,C;Y)& =\sum _{i=1}^{N}logP(Y_i\mid A,B,C)\ge E(A,B,C,\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} l({\mathbf {A}},{\mathbf {B}},{\mathbf {C}}; {\mathbf {Y}})&= \sum _{i=1}^N\log P(Y_{i}\mid {\varvec{A}},{\varvec{B}},{\varvec{C}}) \ge E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi }), \end{aligned}$$\end{document}

with the lower bound $E(A,B,C,\xi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi })$$\end{document} defined as

(11) $\begin{array}{c}E(A,B,C,\xi )=\sum _{i=1}^N{\int }_{{\theta }_i}[\sum _{Z_i}l(A,B,C,{\xi }_i;Y_i,{\theta }_i,Z_i)\times r_i\left(Z_i\right)]\times q_i\left({\theta }_i\right)d{\theta }_i.\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi }) = \sum _{i=1}^N \int _{{\varvec{\theta }}_i}\bigg [\sum _{{\varvec{Z}}_{i}} l({\varvec{A}}, {\varvec{B}},{\varvec{C}}, \varvec{\xi }_i ; Y_i,{\varvec{\theta }}_i,{\varvec{Z}}_{i}) \times r_i({\varvec{Z}}_{i})\bigg ] \times q_i({\varvec{\theta }}_i) d{\varvec{\theta }}_i. \end{aligned}$$\end{document}

Appropriate choices of the variational distributions will lead to a closed form expression of the lower bound in (11). Particularly, following the derivations in Cho et al. (Reference Cho, Wang, Zhang and Xu2021), the above likelihood function implies that an optimal 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$$\end{document} is $q_i\left({\theta }_i\right)\sim N({\theta }_i\mid {\mu }_i,{\Sigma }_i)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ q_i({\varvec{\theta }}_i) \sim N({\varvec{\theta }}_i \mid \mu _i,\Sigma _i )$$\end{document} where the mean and covariance are

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

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

and the variational distributions $r_{\mathrm{ij}}\left(Z_{\mathrm{ij}}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}(Z_{ij})$$\end{document} are $r_{\mathrm{ij}}\left(Z_{\mathrm{ij}}\right)\sim Bernoulli\left(s_{\mathrm{ij}}\right)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}(Z_{ij}) \sim Bernoulli(s_{ij})$$\end{document} , where $s_{\mathrm{ij}}=1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{ij} = 1$$\end{document} if $Y_{\mathrm{ij}}=0$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{ij} = 0$$\end{document} , and otherwise

(14) $\begin{array}{cc}{s}_{\mathrm{ij}}^{-1}=& 1+\frac{c_j}{1-c_j}\frac{1+e^{{\xi }_{i,j}}}{e^{{\xi }_{i,j}}}exp\{-Y_{\mathrm{ij}}({\alpha }_j^{\top }{E}_{q_i}\left[{\theta }_i\right]-b_j)+\\ & \frac{1}{2}(b_j-{\alpha }_j^{\top }{E}_{q_i}\left[{\theta }_i\right]-{\xi }_{i,j})-\eta \left({\xi }_{i,j}\right)\{E_{q_i}\left[{(b_j-{\alpha }_j^{\top }{\theta }_i)}^2\right]-{\xi }_{i,j}^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} s_{ij}^{-1}= & {} 1 + \frac{c_j}{1-c_j} \frac{1+e^{\xi _{i,j}}}{e^{\xi _{i,j}}}\exp \Big \{-Y_{ij}({\varvec{\alpha }}_j^\top E_{q_i}[{\varvec{\theta }}_i]-b_j) + \nonumber \\&\frac{1}{2}(b_j-{\varvec{\alpha }}_j^\top E_{q_i}[{\varvec{\theta }}_i]-{\xi _{i,j}})-\eta ({\xi _{i,j}})\{ E_{q_i}[(b_j-{\varvec{\alpha }}_j^\top {\varvec{\theta }}_i)^2]-\xi _{i,j}^2\} \Big \}. \end{aligned}$$\end{document}

With the above chosen $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$$\end{document} ’s and $r_{\mathrm{ij}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}$$\end{document} ’s, we aim to estimate model parameters $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} , $B$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{B}}$$\end{document} and $C$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{C}}$$\end{document} , together with the introduced local 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} , by maximizing the variational lower bound of the marginal likelihood, $E(A,B,C,\xi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi })$$\end{document} in (11), with the proposed penalties in (6), that is,

(15) $\begin{array}{c}({\hat{A}}_{\lambda },{\hat{B}}_{\lambda },{\hat{C}}_{\lambda },\hat{\xi })={\text{argmax}}_{A,B,C,\xi }E(A,B,C,\xi )-P_{\lambda }\left(A\right)+P\left(B\right)+P\left(C\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} ({\hat{{\mathbf {A}}}}_{\lambda }, {\hat{{\mathbf {B}}}}_{\lambda }, {\hat{{\mathbf {C}}}}_{\lambda },\hat{ \varvec{\xi }}) = \text {argmax}_{{\mathbf {A}},{\mathbf {B}},{\mathbf {C}}, \varvec{\xi }} { E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi }) - P_\lambda ({\mathbf {A}}) + P({\mathbf {B}}) + P({\mathbf {C}})} \end{aligned}$$\end{document}

The corresponding solution $({\hat{A}}_{\lambda },{\hat{B}}_{\lambda },{\hat{C}}_{\lambda })$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\hat{{\mathbf {A}}}}_{\lambda }, {\hat{{\mathbf {B}}}}_{\lambda }, {\hat{{\mathbf {C}}}}_{\lambda })$$\end{document} gives the our GVEM estimators for the penalized likelihood in (6).

To estimate $(A,B,C)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\varvec{A}}, {\mathbf {B}}, {\mathbf {C}})$$\end{document} , we use the coordinate descent algorithm (Friedman, Hastie, Höfling, & Tibshirani, Reference Friedman, Hastie, Höfling and Tibshirani2007; Friedman et al., Reference Friedman, Hastie and Tibshirani2010), which solves the target optimization problem by successively minimizing along each coordinate direction of $(A,B,C)$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\varvec{A}}, {\mathbf {B}}, {\mathbf {C}})$$\end{document} . For each item j, there are one difficulty parameter $b_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} , one guessing parameter $c_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_j$$\end{document} , and K discrimination parameters ${\alpha }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}_j$$\end{document} . The coordinate descent algorithm updates each of the $K+2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K+2$$\end{document} variables according to the following updating rule. (Please see Appendix for a detailed derivation of the updating rule.) Note that the derivation of the below soft-thresholding update rule of $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_{jk}$$\end{document} can be viewed as from the proximal gradient descent algorithm (Beck & Teboulle, Reference Beck and Teboulle2009). Define a function S to be a soft threshold operator such that

(16) $\begin{array}{c}S(\delta ,\lambda )=\text{sign}\left(\delta \right){\left(\right|\delta |-\lambda )}_{+},\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} S(\delta , \lambda ) = \text{ sign }(\delta )(|\delta |-\lambda )_{+}, \end{aligned}$$\end{document}

where for any real number x, sign(x) denotes the sign of x and $x_{+}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_+$$\end{document} denotes $max\{0,x\}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\max \{0,x\}$$\end{document} . The model parameters, ${\alpha }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varvec{\alpha }}_j$$\end{document} ’s, $b_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} and $c_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_j$$\end{document} are updated using Equations (17), (18), and (19), respectively,

(17) $\begin{array}{cc}{\alpha }_{\mathrm{jk}}=& [\sum _{i=1}^N(1-Y_{\mathrm{ij}}+s_{\mathrm{ij}}{Y}_{\mathrm{ij}})(2\eta \left({\xi }_{i,j}\right){[{\Sigma }_i+\left({\mu }_i\right){\left({\mu }_i\right)}^{\top }]}_{k,k}){]}^{-1}\\ & \times S(\sum _{i=1}^N(1-Y_{\mathrm{ij}}+s_{\mathrm{ij}}{Y}_{\mathrm{ij}})\{(Y_{\mathrm{ij}}-\frac{1}{2}){\mu }_{i,k}+2b_j\eta \left({\xi }_{i,j}\right){\mu }_{i,k}\\ & -2\eta \left({\xi }_{i,j}\right)\sum _{l\ne k}{\alpha }_{\mathrm{jl}}{[{\Sigma }_i+\left({\mu }_i\right){\left({\mu }_i\right)}^{\top }]}_{l,k}\},\frac{\lambda }{|{\hat{\alpha }}_{\mathrm{jk}}^{\left(0\right)}{|}^{\gamma }})\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} \alpha _{jk}= & {} \bigg [\sum _{i=1}^N (1-Y_{ij}+s_{ij}Y_{ij})\biggr (2\eta ({\xi _{i,j}})[\Sigma _i+(\mu _i)(\mu _i)^\top ]_{k,k} \biggr )\bigg ]^{-1} \nonumber \\&\times S\biggr (\sum _{i=1}^N(1-Y_{ij}+s_{ij}Y_{ij})\Big \{ (Y_{ij}-\frac{1}{2})\mu _{i,k} + 2b_j\eta ({\xi _{i,j}})\mu _{i,k} \nonumber \\&\qquad -2\eta ({\xi _{i,j}})\sum _{l \ne k}\alpha _{jl}[\Sigma _i+(\mu _i)(\mu _i)^\top ]_{l,k}\Big \} , ~\frac{\lambda }{|{\hat{\alpha }}^{(0)}_{jk}|^\gamma }\biggr ) \end{aligned}$$\end{document}

(18) $\begin{array}{cc}{b}_j=& \frac{{\sum }_{i=1}^N(1-Y_{\mathrm{ij}}+s_{\mathrm{ij}}{Y}_{\mathrm{ij}})[\frac{1}{2}-Y_{\mathrm{ij}}+2\eta \left({\xi }_{i,j}\right){\alpha }_j^{\top }{\mu }_i]+\frac{{\mu }_b}{{\sigma }_b^2}}{2{\sum }_{i=1}^N(1-Y_{\mathrm{ij}}+s_{\mathrm{ij}}{Y}_{\mathrm{ij}})\eta \left({\xi }_{i,j}\right)+\frac{1}{{\sigma }_b^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} b_j= & {} \frac{\sum _{i=1}^N (1-Y_{ij}+s_{ij}Y_{ij})\big [\frac{1}{2}-Y_{ij} + 2\eta (\xi _{i,j}){\varvec{\alpha }}_j^\top \mu _i\big ] +\frac{\mu _b}{\sigma _b^2}}{2\sum _{i=1}^N (1-Y_{ij}+s_{ij}Y_{ij})\eta (\xi _{i,j}) + \frac{1}{\sigma _b^2}}, \end{aligned}$$\end{document}

(19) $\begin{array}{cc}{c}_j=& \frac{{\sum }_{i=1}^N{Y}_{\mathrm{ij}}(1-s_{\mathrm{ij}})+\alpha -1}{N+\alpha +\beta -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} c_j= & {} \frac{\sum _{i=1}^N Y_{ij}(1-s_{ij})+\alpha -1}{N+\alpha +\beta -2}. \end{aligned}$$\end{document}

where ${\hat{\alpha }}_{\mathrm{jk}}^{\left(0\right)}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{\alpha }}^{(0)}_{jk}$$\end{document} is the initial estimator of ${\alpha }_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{jk}$$\end{document} by the GVEM algorithm without including the penalty terms in (15). Additionally, the variational parameter $\xi$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} ’s are updated as

(20) $\begin{array}{cc}{\xi }_{i,j}^2=& b_j^2-2b_j{\alpha }_j^{\top }{\mu }_i+{\alpha }_j^{\top }[{\Sigma }_i+{\mu }_i{\mu }_i^{\top }]{\alpha }_j,\end{array}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \xi _{i,j}^2= & {} b_j^{2} -2b_j{\varvec{\alpha }}_j^{\top } \mu _i + {\varvec{\alpha }}_j^{\top }[\Sigma _i+\mu _i\mu _i^\top ] {\varvec{\alpha }}_j, \end{aligned}$$\end{document}

and the covariance can be updated as

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

To choose the constant sparsity parameter $\lambda$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} , we can apply popular information criteria, such as Akaike Information Criterion (AIC), Bayesian information criterion (BIC) and generalized information criterion (GIC) (Nishii, Reference Nishii1984; Y. Fan & Tang, Reference Fan and Tang2013). We estimate the information criteria by substituting the log-likelihood with the variational lower bound from the GVEM algorithm. The sparsity parameter that minimizes these information criteria will be considered optimal. Our pilot study shows that the GIC method proposed for high-dimensional model selection in Y. Fan and Tang (Reference Fan and Tang2013) performs better than AIC and BIC, and hence GIC is used throughout the study.

The detailed algorithm of the regularized estimation of the loading structure via adaptive Lasso penalization is illustrated in Algorithm 1.

Design

A simulation study was conducted to evaluate the performance of the regularized GVEM algorithm in identifying true item factor loading structure with both M2PL and M3PL models. Three manipulated factors were considered: (1) the number of dimensions was fixed at 3 and 5 (i.e., $K=3,5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K = 3, 5$$\end{document} ); (2) the correlations among factors were fixed at either 0.1, 0.3, or 0.7; and (3) both between-item and within-item multidimensional structures were considered. The sample size was fixed at 2000 (i.e., $N=2000$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N = 2000$$\end{document} ) and 100 replications were run.Footnote 1

For the between-item MIRT model, the test length was 45, with 15 items loaded onto each factor. The true item parameters were selected from the 2013 NAEP item bank (combined national and state assessments) for grade 8. For the within-item MIRT, the true item discrimination parameters were simulated from Unif(0.75, 2), and the difficulty parameters were drawn from the standard normal distribution. Additionally in M3PL, the guessing parameters were fixed at 0.2. The generated item parameters resemble the item parameters in Table 6.1 of (Reckase, Reference Reckase2009) closely. When the dimension was 3, about $60\%$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$60\%$$\end{document} of the items were loaded onto one factor, about $25\%$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$25\%$$\end{document} were loaded onto two factors, and the rest were loaded onto all three factors, whereas for the 5-dimension conditions, about $60\%$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$60\%$$\end{document} , $20\%$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20\%$$\end{document} , $20\%$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20\%$$\end{document} of the items were loaded onto one, two, and three factors, respectively. In all cases, the latent traits $\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} were simulated from $MVN(0,{\Sigma }_{\theta })$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MVN(0,\Sigma _{\varvec{\theta }})$$\end{document} with variance 1, where $r=0.1,0.3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r= 0.1, 0.3$$\end{document} or 0.7.

Six methods were compared in the study, and they are (1) traditional exploratory item factor analysis followed by the CF-Quartimax rotation. This method is denoted as “Rotation” in all results. For this method, during estimation, we did not assume any constraint on the item discrimination parameter but fixed the population covariance matrix to an identity matrix, i.e., ${\Sigma }_{\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{\Sigma }_{\varvec{\theta }}={\varvec{I}}$$\end{document} . The GVEM algorithm was used for model estimation. The final discrimination parameters were transformed to standardized factor loadings, the value of which was compared to 0.3 (Henson & Roberts, Reference Henson and Roberts2006; Costello & Osborne, Reference Costello and Osborne2005). We used ${\psi }_j=U^{-1}{\alpha }_j$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\psi }_j={{\varvec{U}}}^{-1}\varvec{\alpha }_j$$\end{document} , where to obtain the standardized factor loadings. If $|{\psi }_{\mathrm{jk}}|$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\psi _{jk}|$$\end{document} exceeds 0.3, the item is assumed to load on the corresponding factor. This transformation function worked for all simulated conditions except for the within-item structure, $r=0.7,K=3$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7, K=3$$\end{document} , M2PL and M3PL. In these two conditions, we transformed the true discrimination parameters to standardized factor loadings, and found some values were smaller than 0.3. Under these two conditions, we set the cutoff values as 0.75 instead, as the true values were generated from Unif(0.75, 2). Setting a different cutoff will certainly affect the results, and this, to some extent, implies the subjectivity in the traditional EFA rotation method. (2) Exploratory item factor analysis with fixed anchors, and it is denoted as “Fixed Anchors” in all results. For this method, we imposed constraint 1 on the $Q_A$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_A$$\end{document} such that post-hoc rotation is no longer needed. We used the same transformation formula to calculate standardized factor loadings. This method was considered to ensure a direct and fair comparison to the regularization methods. (3) Lasso with constraint 1 and 2; and (4) adaptive Lasso with constraint 1 and 2. For the regularization methods, the tuning parameter $\lambda$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} was chosen by GIC. The GIC was computed as follows:

$\begin{array}{c}GIC=log(log(N\left)\right)\times log\left(N\right)\times k-2\times E(A,B,C,\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} GIC = \log (\log (N))\times \log (N) \times k -2\times E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi }), \end{aligned}$$\end{document}

where N refers to the sample size, k refers to the number of parameters estimated by the model, $E(A,B,C,\xi )$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E({\varvec{A}},{\varvec{B}},{\varvec{C}}, \varvec{\xi })$$\end{document} refers to the lower bound.

In additional to the two constraints for the model ability, we truncated ${\hat{\alpha }}_{\mathrm{jk}}$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{{\varvec{\alpha }}}}_{jk}$$\end{document} to 0 if $|{\hat{\alpha }}_{\mathrm{jk}}|<0.001$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{\hat{{\varvec{\alpha }}}}_{jk}| < 0.001$$\end{document} . As to $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} in adaptive Lasso, Zou (2006) recommended three values, 0.5, 1, and 2. A few pilot trials were conducted to decide on the optimal $\gamma$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , and $\gamma =2$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =2$$\end{document} was used for all conditions except a few conditions in which case $\gamma =1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =1$$\end{document} was used. These conditions are within-item M2PL, $r=0.7,k=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7, k=5$$\end{document} , constraint 1 and 2, as well as between-item M2PL, $r=0.7,k=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7, k=5$$\end{document} , constraint 2 only.

As the main objective of this section is to estimate relationship between test items and latent traits, we used the correct estimation rate of A matrix (eq. (22)). It measures how well the sparsity of the A matrix is estimated by the regularized estimation. Notice that we only calculated correct rate for entries excluding the first K by K sub-matrix since we fixed this part to have identity matrix as a zero structure to ensure identifiability.

(22) $\begin{array}{c}CR=\frac{1}{K\times J}\sum _{1\le j<J,1\le k\le K}I({\hat{Q}}_{\mathrm{jk}}=Q_{\mathrm{jk}}^{\mathrm{true}})\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} CR = \frac{1}{K\times J} \sum _{1 \le j < J, 1 \le k \le K} I({{\hat{Q}}}_{jk} = Q^{true}_{jk}) \end{aligned}$$\end{document}

We also compared the performance of Lasso and adaptive Lasso penalization using two measures: sensitivity and specificity. In our context, sensitivity is the probability of correctly identifying nonzero entries among true nonzero entries. Specificity is the probability of correctly identifying zero entries among true zero entries. In other words, sensitivity measures the true negative rate, while specificity illustrates the true positive rate. Naturally, a test with both high sensitivity and high specificity is desired, although there is always a trade-off.

Other criteria include the average relative bias and root mean squared error (RMSE). The parameter recovery for ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Sigma _{\varvec{\theta }}$$\end{document} is calculated by taking differences between each freely estimated entries of the true ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Sigma _{\varvec{\theta }}$$\end{document} and estimated ${\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{\Sigma }}_{\varvec{\theta }}$$\end{document} . Relative bias and RMSE were obtained for each nonzero model parameter across all items within a condition first and then averaged over 100 replications.

Simulation Results

In this section, we first present the simulation results under various settings in M2PL and M3PL with boxplots to show the distribution of correct estimation rates, sensitivities, and specificities. Among the three information criteria, GIC showed the best performance at selecting the optimal result as it favors the models that penalizes more on the number of parameters; thus, we present the simulation results with GIC selection criteria in figures in this section.

Figures 1 and 2 show the recovery of item factor loading structure in terms of correct rates, sensitivity and specificity under M2PL and M3PL, respectively. All six methods were presented in the same order under each manipulated condition in Figures 1, 2, 3, 4, 5, 6. For M2PL, the adaptive Lasso method is consistently the best-performing method under all conditions, except when $r=0.1,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.1, K=5$$\end{document} , item structure is within-item M2PL, and when $r=0.3,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.3, K=5$$\end{document} , item structure is between-item M2PL. Under these two conditions, EFA rotation method performs slightly better than the adaptive Lasso method. The EFA with fixed anchors and Lasso regularization methods, on the other hand, performs a lot worse. When $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 within-item M2PL is used, EFA rotation method performs considerably worse than EFA with fixed anchors and adaptive Lasso methods. Between the two constraint settings, constraint 2 yields more free parameters and hence it is harder to handle than constraint 1. Therefore, it is not surprising that adaptive Lasso with constraint 1 performs slightly better than with constraint 2 in more challenging scenarios (i.e., higher correlation, larger K, and within-item multidimensionality), whereas the difference between the two types is almost negligible in simpler scenarios.

When M3PL is the data generating model, the recovery of item factor loading structure is generally worse than that from M2PL, with a decrement of correct rate, sensitivity, and specificity in the range of 5% to 20%. The general trend of the manipulated factors on the results stay the same as compared to M2PL. That is, increasing factor correlation or allowing item cross-loadings makes the recovery of factor structure harder, although adaptive Lasso still performs the best among the six methods in all conditions except when $r=0.7$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7$$\end{document} , $K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=5$$\end{document} and test exhibits between-item multidimensional structure. In this case, EFA rotation method tends to excel.

Figure. 1

Correct estimation rates of item factor loading structure under M2PL.

Figure. 2

Correct estimation rates of item factor loading structure under M3PL.

Figure 3 presents the relative bias of model parameters under M2PL. When the test has 5 latent factors and $r=0.1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.1$$\end{document} or 0.3, although relative bias vary slightly differently for different parameters, the results from the six methods are almost indistinguishable. In a between-item structure with 3 factors, the relative bias for b has more variability across replications. It is because the true parameters of some items are close to 0. The relative bias vary more for the within-item condition in general. In a within-item structure, the two regularization methods appear to produce less bias than the EFA rotation method especially for ${\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} . Under $K=3,r=0.1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3, r=0.1$$\end{document} or 0.7, within-item M2PL conditions, the relative bias values for ${\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} estimated by EFA rotation fall outside of the range. Figure 4 presents the RMSE of model parameters under M2PL. Again, all six methods produce comparable RMSE when $r=0.1$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.1$$\end{document} under between-item condition. When the factor correlation increases, the EFA rotation method generates larger RMSE for $\alpha$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} and ${\Sigma }_{\theta }$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Sigma _\theta $$\end{document} . The same trend holds under the within-item M2PL conditions, although adaptive Lasso method seems to generate large RMSE for some conditions. Under the most difficult condition of $r=0.7,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7, K=5$$\end{document} , one can see a lot of variability of RMSE across replications. In this case, the two Lasso methods seem to produce smaller median RMSE for majority of the parameters than EFA rotation method and EFA with fixed anchors method. The better performance of Lasso methods compared to adaptive Lasso may be because of two reasons: (1) we computed bias and RMSE only on those parameters whose true values were nonzero. Hence, even if the Lasso method fails to shrink some true zero loadings to zero, they will not count toward bias or RMSE. (2) Initial values play an important role in adaptive Lasso to determine an adaptive penalty weight. We used the results from EFA rotation methods as initial values, and other better initial values could be explored in the future, such as the SVD method in H. Zhang et al. (Reference Zhang, Chen and Li2020).

Figure. 3

Relative bias of model parameter estimates under M2PL.

Figure. 4

RMSE of model parameter estimates under M2PL.

Figure. 5

Relative bias of model parameter estimates under M3PL.

Figure. 6

RMSE of model parameter estimates under M3PL.

Figures 5 and 6 show the relative bias and RMSE of model parameters under M3PL. The inclusion of the guessing parameter, unsurprisingly, makes the model parameter recovery much harder, as shown in larger bias and RMSE as well as more variability across replications. The overall pattern observed from M2PL results continued to hold. That is, increasing factor correlation and using within-item factor structure not only increase relative bias and RMSE, but also yield more instability across replications. The EFA rotation method produces the largest average absolute bias and mean RMSE in almost all conditions, followed by EFA with fixed anchors method, although results from regularization methods seem to have more variability when the factor correlation is high.

In summary, the adaptive Lasso method outperforms EFA rotation method under almost all conditions regarding item factor loading structure recovery. There are only 3 exceptions: when $r=0.1,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.1, K=5$$\end{document} , item structure is within-item M2PL, when $r=0.3,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.3, K=5$$\end{document} , item structure is between-item M2PL, and $r=0.7,K=5$ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=0.7, K=5$$\end{document} , item structure is between-item M3PL. In these three conditions, EFA rotation method performs better than the adaptive Lasso method by a small margin. Under some simple scenarios (i.e., low-correlation or medium-correlation and $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} , item structure is between-item M2PL), there is no appreciable difference between the EFA rotation method and the adaptive Lasso method with either type of constraints. As for item parameter recovery, the adaptive Lasso method outperforms EFA rotation method for all of the high-correlation scenarios in M2PL. For small-correlation, between-item M2PL conditions, the results of adaptive Lasso and EFA rotation method appear to be indistinguishable. In M3PL, the adaptive Lasso method produces more accurate results compared to EFA rotation method under all conditions. Only under between-item M3PL conditions, EFA rotation generates smaller RMSE values and relative bias with less variability for ${\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} , but it produces larger RMSE and relative bias for other parameters.