2.1 Log-linear cognitive diagnostic model

In this study, we focused on the LCDM. This is because it is a general model that contains a large number of models that have been previously discussed, such as DINA, DINO, rRUM, and LLM (Henson et al., Reference Henson, Templin and Willse2009). More importantly, the LCDM can provide a parameterization that not only enables it to characterize the differences between the various models but also offers support for more complex data structures (Henson et al., Reference Henson, Templin and Willse2009). In fact, any possible set of constraints for the saturated form LCDM can be used to define a model that fits the item response in the framework of cognitive theory. Moreover, a better understanding of the relationships between compensatory models and non-compensatory models can be described in the general parametric form. After this, a brief introduction to the LCDM will be given.

First, we define several indices that will be important throughout this paper. Each examinee is denoted by $i~(i=1,\cdots ,N)$ , each item by $j~(j=1,\cdots ,J)$ , each attribute by $k~(k=1,\cdots ,K)$ , and the latent class corresponding to an attribute profile is denoted by $l~(l=1,\cdots ,L)$ . We consider the latent attribute $\alpha _{ik}$ to be a binary variable, where the absence or presence of the corresponding attribute is represented by the values 0 and 1, respectively. $\boldsymbol {\alpha }_{i}=(\alpha _{i1},\cdots ,\alpha _{ik},\cdots ,\alpha _{iK})^{\text {T}}$ is a vector of K-dimensional latent attribute profiles for the ith examinee. In light of the categorical nature of the latent classes, $\boldsymbol {\alpha }_{i}$ belongs to one of $L=2^K$ latent attribute profiles. Defining $\boldsymbol {\widetilde {\alpha }}_{l}=(\widetilde {\alpha }_{l1},\cdots ,\widetilde {\alpha }_{lk},\cdots ,\widetilde {\alpha }_{lK})^{\text {T}}$ as the attribute profile for examinees of class l, where $\widetilde {\alpha }_{lk}$ is 1 if the examinees of class l acquire skill k and 0 otherwise, will be useful in the following.

$\widetilde {\textbf {A}}=(\widetilde {\boldsymbol {\alpha }}_{1},\cdots ,\widetilde {\boldsymbol {\alpha }}_{l},\cdots ,\widetilde {\boldsymbol {\alpha }}_{L})^{\text {T}}$ denotes a matrix of $L \times K$ dimensions containing all the attribute profiles. The $\text {Q}$ -matrix(Tatsuoka, Reference Tatsuoka1983) is a $J \times K $ matrix used to describe the relationship between attributes and items, where $\boldsymbol {q}_j^{\text {T}}=(q_{j1},\cdots ,q_{jk},\cdots ,q_{jK})$ , and $q_{ik}\in \{0,1\}$ is a vector of the jth row of the $\text {Q}$ -matrix; that is, $\textbf {Q}=(\boldsymbol {q}_1,\cdots ,\boldsymbol {q}_j,\cdots ,\boldsymbol {q}_J)^{\text {T}}$ : $q_{ik}=1$ if the attribute k is required by item j, and $q_{ik}=0$ otherwise. Next, a binary latent indicator variable $\boldsymbol {z}_i=[z_{i1},\cdots ,z_{il},\cdots ,z_{iL}]^{\text {T}}$ is introduced, which satisfies $\sum _{l=1}^{L}z_{il}=1$ , where $z_{il}=1$ denotes the ith examinee belonging to the lth attribute profile (i.e., $\boldsymbol {\alpha }_{i}=\widetilde {\boldsymbol {\alpha }}_{l}$ ). Let $x_{ij}$ be the observed item response for the ith examinee to the jth item: $x_{ij}=1$ if the ith examinee gives the correct answer for the jth item, and it is 0 otherwise. The corresponding item response matrix for all examinees answering all items is $\textbf {X}=(\boldsymbol {x}_{1},\cdots ,\boldsymbol {x}_{i},\cdots ,\boldsymbol {x}_{N})^{\text {T}}$ , where $\boldsymbol {x}_{i}=(x_{i1},\cdots ,x_{ij},\cdots ,x_{iJ})^{\text {T}}$ , $i=1,\cdots ,N$ . Then, the probability of a correct response for the LCDM can be expressed as

(1) $$ \begin{align} P(x_{ij}=1|\boldsymbol{\alpha}_{i}=\widetilde{\boldsymbol{\alpha}}_{l},\eta_j,\boldsymbol{\lambda}_j,\boldsymbol{q}_j)=\frac{\exp(\boldsymbol{\lambda}_{j}^{\text{T}}\textbf{h}(\widetilde{\boldsymbol{\alpha}}_{l},\boldsymbol{q}_{j})+\eta_{j})}{1+\exp(\boldsymbol{\lambda}_{j}^{\text{T}}\boldsymbol{h}(\widetilde{\boldsymbol{\alpha}}_{l},\boldsymbol{q}_{j})+\eta_{j})}, \end{align} $$

where $\eta _{j}$ is the intercept parameter, and $\exp (\eta _{j})/(1+\exp (\eta _{j}))$ indicates the probability that an examinee answers correctly on item j if he or she does not master any of the attributes examined on that item. $\boldsymbol {\lambda }_j=(\lambda _{j1},\cdots ,\lambda _{jK},\lambda _{j12},\cdots ,\lambda _{j12\cdots (K-1)K})^{\text {T}}$ is the slope parameter vector, which is composed of a $D \times 1$ vector, where $D=2^K-1$ . $\textbf {h}(\widetilde {\boldsymbol {\alpha }}_{l},\boldsymbol {q}_{j})$ represents a set of linear combinations of $\widetilde {\boldsymbol {\alpha }}_{l}$ and $\boldsymbol {q}_{j}$ :

(2) $$ \begin{align} \boldsymbol{\lambda}_{j}^{\text{T}}\textbf{h}(\widetilde{\boldsymbol{\alpha}}_{l},\boldsymbol{q}_{j})=\sum_{k=1}^{K}\lambda_{jk}\widetilde{\alpha}_{lk}q_{jk}+\sum_{k=1}^{K}\sum_{k^{\prime}>k}\lambda_{jkk^{\prime}}\widetilde{\alpha}_{lk}\widetilde{\alpha}_{ik^{\prime}}q_{jk}q_{jk^{\prime}}+\cdots+\lambda_{j12\cdots(K-1)K}\prod_{k=11}^{K}\widetilde{\alpha}_{lk}q_{jk}. \end{align} $$

Combining the latent variable $\boldsymbol {z}_i$ and Eq. (2), the LCDM can be rewritten as

(3) $$ \begin{align} P(x_{ij}=1|\widetilde{\textbf{A}},\boldsymbol{z}_i,\eta_j,\boldsymbol{\lambda}_j,\boldsymbol{q}_j)=\prod_{l=1}^LP(x_{ij}=1|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\eta_j,\boldsymbol{\lambda}_j,\boldsymbol{q}_j)^{z_{il}}. \end{align} $$

2.2 DINA model

The DINA model, as a special case of the LCDM, has a relatively straightforward structure and widespread adoption in cognitive diagnostic assessments; specialized software packages are also available for a number of estimation techniques grounded in the model. Therefore, we provide a short overview of the traditional DINA model and its interconversion with the LCDM. Two-item parameters have been introduced in the traditional DINA models for each item j: $s_j$ is the slipping parameter and $g_j$ is the guessing parameter, and the probability of a correct response can be written as

(4) $$ \begin{align} \begin{aligned} P(x_{ij}=1|\boldsymbol{\alpha}_i=\widetilde{\boldsymbol{\alpha}}_{l},g_j,s_j)&=g_j^{1-\gamma_{lj}}(1-s_j)^{\gamma_{lj}},\\ \gamma_{lj}&=\prod_{k=1}^{K}\alpha_{lk}^{q_{jk}}, \end{aligned} \end{align} $$

where $\gamma _{ij}$ is the ideal response pattern. $\gamma _{lj}=1$ indicates that examinee i possesses all the required attributes for item j; otherwise, $\gamma _{lj}=0$ . The parameters $s_j$ and $g_j$ can be formally defined by

(5) $$ \begin{align} \begin{aligned} s_j&=P(x_{ij}=0|\gamma_{lj}=1),\\ g_j&=P(x_{ij}=1|\gamma_{lj}=0). \end{aligned} \end{align} $$

Since the estimation approach presented in this work is based on the LCDM, we must first convert the DINA model to LCDM format. Our next topic is the connection between the DINA model and the LCDM and how they may be converted back and forth.

Let $\widetilde {\boldsymbol {K}}_j=\{k: \text {attribute}~ k ~\text {is measured by item}~ j\}$ denote an indicator set of attributes investigated by item j and $K^*_j$ denote the number of investigated attributes. Then, the DINA model can be rewritten in the form:

(6) $$ \begin{align} P(x_{ij}=1|\boldsymbol{\alpha}_i=\widetilde{\boldsymbol{\alpha}}_{l},g_j,s_j)=\frac{\exp\left(\eta_{j}+\lambda_{j\tilde{K}_{j1}\tilde{K}_{j2}\cdots \tilde{K}_{jK_j^*}}\prod\limits_{k^*\in \tilde{K}_j}^{}\alpha_{lk^*}q_{jk^*}\right)}{1+\exp\left(\eta_{j}+\lambda_{j\tilde{K}_{j1}\tilde{K}_{j2}\cdots \tilde{K}_{jK_j^*}}\prod\limits_{k^*\in \tilde{K}_j}^{}\alpha_{lk^*}q_{jk^*}\right)}, \end{align} $$

where

(7) $$ \begin{align} \begin{aligned} \eta_j&=-\log \left(\frac{g_j}{1-g_j}\right),\\ \lambda_{j\tilde{K}_{j1}\tilde{K}_{j2}\cdots \tilde{K}_{jK_j^*}}&=-\eta_{j}+\log \left(\frac{1-s_j}{s_j}\right). \end{aligned} \end{align} $$

For simplicity, we denote $\lambda _{j\tilde {K}_{j1}\tilde {K}_{j2}\cdots \tilde {K}_{jK_j^*}}$ as $\lambda _j$ and the DINA model is equivalent to the following form:

(8) $$ \begin{align} P(x_{ij}=1|\boldsymbol{\alpha}_i=\widetilde{\boldsymbol{\alpha}}_{l},g_j,s_j)=\frac{\exp\left(\eta_{j}+\lambda_{j}\prod\limits_{k=1}^{K}\alpha_{lk}^{q_{jk}}\right)}{1+\exp\left(\eta_{j}+\lambda_{j}\prod\limits_{k=1}^{K}\alpha_{lk}^{q_{jk}}\right)}. \end{align} $$

While this study focuses mostly on the LCDM, various variants of the LCDM, such as the DINO model, LLM, and saturated LCDM, are also discussed. We will therefore not go into great depth here; instead, the reader should refer to the Supplementary Material for the necessary information.

3.1 Variational Bayesian EM algorithm

Since it is straightforward to convert an approximate conditional posterior distribution problem into an optimization problem using VI methods, these techniques see extensive application in inferring Bayesian models in the area of machine learning (Beal, Reference Beal2003; Bishop, Reference Bishop2006; Jordan et al., Reference Jordan, Ghahramani, Jaakkola and Saul1999). Next, we briefly outline the implementation process of the variational Bayesian EM (VBEM) algorithm (Beal, Reference Beal2003). Assume that the observed dataset $\boldsymbol {y}=(\boldsymbol {y}_1,\cdots ,\boldsymbol {y}_i,\cdots ,\boldsymbol {y}_N)$ is produced by model $\mathcal {M}$ , where model $\mathcal {M}$ consists of the latent variables $\boldsymbol {\zeta }=(\boldsymbol {\zeta }_1,\cdots ,\boldsymbol {\zeta }_{i},\cdots ,\boldsymbol {\zeta }_{N}$ ) and model parameters $\boldsymbol {\theta }$ . Next, we specify a variational density family $\mathcal {Q}$ over the unknown variables $\boldsymbol {\zeta }$ and $\boldsymbol {\theta }$ . The purpose of this is to establish the optimal approximation $q(\boldsymbol {\zeta },\boldsymbol {\theta })\in \mathcal {Q}$ to their posterior distribution using this specified variational density (i.e., $q(\boldsymbol {\zeta },\boldsymbol {\theta })\rightsquigarrow p(\boldsymbol {\zeta },\boldsymbol {\theta }|\boldsymbol {y})$ ). Next, we introduce the concept of the evidence lower bound (ELBO), which is critical for determining the optimal $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ . Let $p(\boldsymbol {y}|\mathcal {M})$ be a marginal density of the model $\mathcal {M}$ ; the ELBO can then be represented as a lower bound of the logarithm marginal density $\log p(\boldsymbol {y}|\mathcal {M})$ :

(9) $$ \begin{align} \begin{aligned} \log p(\boldsymbol{y}|\mathcal{M})&=\log \int p(\boldsymbol{y},\boldsymbol{\zeta},\boldsymbol{\theta})\,d\boldsymbol{\zeta}d\boldsymbol{\theta}\\ &=\log \int q(\boldsymbol{\zeta},\boldsymbol{\theta}) \frac{p(\boldsymbol{y},\boldsymbol{\zeta},\boldsymbol{\theta})}{q(\boldsymbol{\zeta},\boldsymbol{\theta})}\, d\boldsymbol{\zeta}d\boldsymbol{\theta}\\ &\geq \int q(\boldsymbol{\zeta},\boldsymbol{\theta})\log \frac{p(\boldsymbol{y},\boldsymbol{\zeta},\boldsymbol{\theta})}{q(\boldsymbol{\zeta},\boldsymbol{\theta})}\, d\boldsymbol{\zeta}d\boldsymbol{\theta}~~~~~\textbf{Jensen's inequality}\\ &=\mathrm{E}_{q(\boldsymbol{\zeta},\boldsymbol{\theta})} [\log p(\boldsymbol{y},\boldsymbol{\zeta},\boldsymbol{\theta})-\log q(\boldsymbol{\zeta},\boldsymbol{\theta})]\\ &\triangleq \underline{\mathcal{L}}\left(q(\boldsymbol{\zeta},\boldsymbol{\theta})\right) ,\\ \end{aligned} \end{align} $$

where $\underline {\mathcal {L}}\left (q(\boldsymbol {\zeta },\boldsymbol {\theta })\right )$ is denoted as the ELBO, which is a function of the free distribution $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ . We need to maximize $\underline {\mathcal {L}}\left (q(\boldsymbol {\zeta },\boldsymbol {\theta })\right )$ with respect to $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ so that it tends more closely to $ \log p(\boldsymbol {y}|\mathcal {M})$ . Blei et al. (Reference Blei, Kucukelbir and McAuliffe2017) presented a formula connecting $ \log p(\boldsymbol {y}|\mathcal {M})$ with the ELBO and the Kullback–Leibler (KL) divergence:

(10) $$ \begin{align} \log p(\boldsymbol{y}|\mathcal{M})=\underline{\mathcal{L}}\left(q(\boldsymbol{\zeta},\boldsymbol{\theta})\right)+\mathrm{KL}(q(\boldsymbol{\zeta},\boldsymbol{\theta}) \Vert p(\boldsymbol{\zeta},\boldsymbol{\theta}|\boldsymbol{y})). \end{align} $$

Since $\log p(\boldsymbol {y}|\mathcal {M})$ is a constant with respect to $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ , maximizing the ELBO is actually equivalent to minimizing the KL distance. Specifically, the optimal $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ we obtained in the variational density family $\mathcal {Q}$ is the density that minimizes the KL divergence between the posterior distribution $p(\boldsymbol {\zeta },\boldsymbol {\theta }|\boldsymbol {y})$ and itself. To further simplify the variational density $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ , we assume that it satisfies mean-field theory. Mean-field theory has been widely used in variational Bayesian inference (Beal, Reference Beal2003; Blei et al., Reference Blei, Kucukelbir and McAuliffe2017; Jordan et al., Reference Jordan, Ghahramani, Jaakkola and Saul1999; Wand et al., Reference Wand, Ormerod, Padoan and Frühwirth2011). In the mean-field theory, latent variables are mutually independent and each is governed by a separate factor in the variational density, allowing the variational density $q(\boldsymbol {\zeta },\boldsymbol {\theta })$ to be decomposed into $q(\boldsymbol {\zeta })q(\boldsymbol {\theta })$ . An iterative optimization procedure is implemented by seeking to maximize the mean-field variational density of a parameter of interest while fixing the others. The VB algorithm can be divided into the following two steps:

(11) $$ \begin{align} \begin{aligned} \textbf{VBE~step:} &\qquad\qquad q^{\mathrm{new}}(\boldsymbol{\zeta}_i)=\frac{1}{\mathcal{Z}_{\boldsymbol{\zeta}_i}}\exp \left[ \int q^{\mathrm{old}}(\boldsymbol{\theta})\log p(\boldsymbol{\zeta}_i,\boldsymbol{y}_i|\boldsymbol{\theta})\,d\boldsymbol{\theta}\right]\\ &\qquad \qquad \qquad \qquad \propto\exp \left\{\mathrm{E}_{q^{\mathrm{old}}(\boldsymbol{\theta})}[\log p(\boldsymbol{y}_i,\boldsymbol{\zeta},\boldsymbol{\theta})]\right\}, for~\forall i, \\ \textbf{VBM~step:}&\qquad\qquad q^{\mathrm{new}}(\boldsymbol{\theta})=\frac{1}{\mathcal{Z}_{\boldsymbol{\theta}}}p(\boldsymbol{\theta})\exp \left[ \int q^{\mathrm{new}}(\boldsymbol{\zeta}) \log p(\boldsymbol{\zeta},\boldsymbol{y}|\boldsymbol{\theta})\,d\boldsymbol{\zeta}\right]\\ &\qquad \qquad \qquad \qquad \propto \exp \left\{\mathrm{E}_{q^{\mathrm{new}}(\boldsymbol{\zeta})}[\log p(\boldsymbol{y},\boldsymbol{\zeta},\boldsymbol{\theta})]\right\}, \end{aligned} \end{align} $$

where $\mathcal {Z}_{\boldsymbol {\zeta }_i}$ and $\mathcal {Z}_{\boldsymbol {\theta }}$ are the normalizing constants. To sum up, the variational density for the latent variable is updated in the VBE step, while the variational density for the model parameters is updated in the VBM step. Therefore, the prerequisite to be able to implement the VBEM algorithm is that the posterior distribution of all parameters, either latent variables or model parameters, should have a closed form. The VEM algorithm proposed by Yamaguchi and Okada (Reference Yamaguchi and Okada2020a, Reference Yamaguchi and Okada2020b) in educational psychometric research is essentially identical to the VBEM algorithm provided by Beal (Reference Beal2003), with the only differences being in nomenclature.

3.2 Variational methods in Bayesian logistic regression

As mentioned above, implementing the VBEM algorithm requires a closed form for the posterior distributions of each parameter. Therefore, the VBEM algorithm cannot be directly applied to the LCDM based on the logit link function. To overcome this challenge, we adopt Jaakkola and Jordan’s (Reference Jaakkola and Jordan2000) variational Bayesian method for logistic regression models to estimate the more complex LCDM in the cognitive diagnostic framework. Specifically, their method uses a Taylor expansion on the logistic function to obtain a tight lower bound, facilitating parameter representation in a Gaussian distribution form that is easily implementable for VI. Next, we will provide the mathematical expression that Jaakkola and Jordan (Reference Jaakkola and Jordan2000) used for performing the first-order Taylor expansion and the specific derivation of the tight lower bound.

Consider the logistic function $\sigma (\omega )=1/\left (1+\exp (-\omega )\right )$ . The corresponding log logistic function can be derived as

(12) $$ \begin{align} \log \sigma(\omega)=-\log \left(1+\exp(-\omega)\right)=\frac{\omega}{2}-\log \left(\exp\left(\frac{\omega}{2}\right)+\exp\left(-\frac{\omega}{2}\right)\right). \end{align} $$

Denote that

$$ \begin{align*}f(\omega)=-\log \left(\exp\left(\frac{\omega}{2}\right)+\exp\left(-\frac{\omega}{2}\right)\right). \end{align*} $$

By calculating the second derivative, we can determine that $f(\omega )$ is a convex function about the variable $\omega ^2$ . Therefore, any tangent line of $f(\omega )$ can serve as its lower bound, as it will always be less than or equal to $f(\omega )$ . A tight lower bound function for $f(\omega )$ can be obtained by executing a first-order Taylor expansion on the function $f(\omega )$ in terms of the variable $\omega ^2$ at the point $\xi ^2$ ,

(13) $$ \begin{align} f(\omega) \geq f(\xi)+\frac{\partial f(\xi)}{(\partial\xi^2)}(\omega^2-\xi^2)=f(\xi)-\frac{1}{2\xi}(\sigma(\xi)-\frac{1}{2})(\omega^2-\xi^2). \end{align} $$

According to Eqs.(12) and (13), we can derive a tight lower bound of $\sigma (\omega )$ with the specific form as

(14) $$ \begin{align} \sigma(\omega)\geq\sigma(\xi) \exp\left(\frac{(\omega-\xi)}{2}-\tau(\xi)(\omega^2-\xi^2)\right),\quad \tau(\xi)=\frac{1}{2\xi}\left(\sigma(\xi)-\frac{1}{2}\right), \end{align} $$

which results in a quadratic form on $\omega $ .

Regarding the LCDM, which also employs a logistic form, $\omega $ represents a set of linear combinations. These combinations involve unknown item parameters, an individual’s latent attribute vector, and the known $\text {Q}$ -matrix within the LCDM framework (for further details, please refer to Eqs. (1) and (2)). Based on Eq. (14), we can derive a quadratic form for the item parameter. Consequently, the VI algorithm can be implemented in the LCDM by deriving a specific posterior distribution for item parameters using the Gaussian prior distribution, which serves as the conjugate prior for these parameters. In the next subsection, we will focus on elucidating the process of deriving the tight lower bound in the LCDM using Eq. (14).

3.3 Tight lower bound for the LCDM

In this section, the goal is to derive the tight lower bound for LCDM as outlined above. We first conducted a transformation on the item response data in the LCDM to make it easier to acquire the tight lower bound term of the likelihood function before providing the implementation of our VBEM-M algorithm. The item response data $x_{ij}=\{0,1\}$ is transformed into $y_{ij}=\{-1,1\}$ with the help of the equation $y_{ij}=2x_{ij}-1$ . Let $\boldsymbol {\lambda }^*_j=(\eta _j,\boldsymbol {\lambda }_j^{\text {T}})^{\text {T}}$ and $\boldsymbol {h}_{jl}^*=(1,h(\boldsymbol {\tilde {\alpha }}_{l},\boldsymbol {q}_j))$ ; the item response probability of $y_{ij}$ is then given by

(15) $$ \begin{align} \begin{aligned} P(y_{ij}=1|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\boldsymbol{\lambda}^*_j,\boldsymbol{q}_j)=P(x_{ij}=1|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\boldsymbol{\lambda}^*_j,\boldsymbol{q}_j) =\frac{1}{1+\exp(-\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*)},\\ P(y_{ij}=-1|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\boldsymbol{\lambda}_j^*,\boldsymbol{q}_j)=P(x_{ij}=0|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\boldsymbol{\lambda}_j^*,\boldsymbol{q}_j)=\frac{1}{1+\exp(\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}^*_{jl})}. \end{aligned} \end{align} $$

Recalling the logistic function form, the item response probability of $y_{ij}$ can then be rewritten as follows:

(16) $$ \begin{align} \begin{aligned} &p(y_{ij}|\boldsymbol{\alpha}_i=\tilde{\boldsymbol{\alpha}_l},\boldsymbol{\lambda}^*_j,\boldsymbol{q}_j)=\sigma(y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}^*_{jl}).\\ \end{aligned} \end{align} $$

Therefore, the likelihood based on the introduced latent variable z can be represented by

(17) $$ \begin{align} p(\textbf{Y}|\boldsymbol{z},\widetilde{\textbf{A}},\boldsymbol{\lambda}^*,\textbf{Q})=\prod_{i=1}^N\prod_{j=1}^J \prod_{l=1}^L\sigma(y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}^*_{jl})^{z_{il}}. \end{align} $$

According to Eq. (14), the tight lower bound function for $\sigma (y_{ij}\boldsymbol {\lambda }_j^{*\text {T}}\boldsymbol {h}^*_{jl})$ is determined by performing a first-order Taylor expansion with respect to the variable $(y_{ij}\boldsymbol {\lambda }_j^{*\text {T}}\boldsymbol {h}^*_{jl})^2$ at the point $\xi _{ijl}^2$ . Therefore, a tight lower bound of the likelihood for the LCDM can be derived by:

(18) $$ \begin{align} \begin{aligned} &p(\textbf{Y}|\boldsymbol{z},\widetilde{\textbf{A}},\boldsymbol{\lambda}^*,\textbf{Q})=\prod_{i=1}^N\prod_{j=1}^J \prod_{l=1}^L\sigma(y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}^*_{jl})^{z_{il}}\\ &\geq \prod_{i=1}^N\prod_{j=1}^J \prod_{l=1}^L \left\{\sigma(\xi_{ijl})\exp\left(\frac{y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*-\xi_{ijl}}{2}-\tau(\xi_{ijl})(\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*\boldsymbol{h}_{jl}^{*\text{T}}\boldsymbol{\lambda}_j^*-\xi_{ijl}^{2})\right)\right\}^{z_{il}}\\ &\triangleq \underline{p(\textbf{Y}|\boldsymbol{z},\widetilde{\textbf{A}},\boldsymbol{\lambda}^*,\textbf{Q})}. \end{aligned} \end{align} $$

Given that the tight lower bound for the likelihood function is of exponential form, using the multivariate normal distribution as a conjugate prior distribution for $\boldsymbol {\lambda }^*$ will yield a closed-form posterior distribution. Due to these considerations, in the subsequent computations, we implement the VBEM-M algorithm using the tight lower bound of the likelihood function rather than the original likelihood function. Moreover, it is important to highlight that a new local parameter, $\xi _{ijl}$ , has been introduced at this stage. Determining the optimal value for $\xi _{ijl}$ is an essential part of our analysis. In this paper, we implement a maximization process to ascertain the most suitable value for $\xi _{ijl}$ . The detailed methodology behind this process will be elaborated in the following subsection.

3.4 Fully Bayesian representation of the joint posterior distribution

In the fully Bayesian framework, statistical inference relies on the selection of the prior distribution. The posterior distribution can be derived by combining the prior distribution (prior information) with the likelihood function (sample information). Prior distributions from the following Bayesian hierarchical structures will be considered in this study:

(19) $$ \begin{align} \begin{aligned} & y_{ij} \sim p( y_{ij}|\boldsymbol{z}_i,\widetilde{\textbf{A}},\boldsymbol{\lambda}^*,\textbf{Q}),~p\left(\boldsymbol{z}_i | \boldsymbol{\pi}\right) = \prod_{l=1}^L\pi_l^{z_{il}}, 0 \leq \pi_l\leq 1, \sum_{l=1}^L \pi_l=1,\\ & p(\boldsymbol{\pi})=p\left(\pi_1, \ldots, \pi_L\right) ~=~ \text {Dirichlet }\left(\boldsymbol{\delta}_0\right), \boldsymbol{\delta}_0=\left(\delta_{01}, \ldots, \delta_{0L}\right),\\ & p(\boldsymbol{\lambda}_{j}^*)~=~ \text {MVN}\left(\boldsymbol{\lambda}_0^*,\textbf{I}_{D+1}\right), \boldsymbol{\lambda}_0^*=(\eta_0,\boldsymbol{\lambda}_0), \\ &\boldsymbol{\lambda}_0 ~=~(\underbrace{\lambda_{0,1},\cdots,\lambda_{0,K}}_{\textbf{main effect terms}},\underbrace{\lambda_{0,K+1},\cdots,\lambda_{0,D}}_{\textbf{interaction terms}}) ~=~(\underbrace{\lambda_{0,main},\cdots,\lambda_{0,main}}_{\textbf{main effect terms}},\underbrace{\lambda_{0,inter},\cdots,\lambda_{0,inter}}_{\textbf{interaction terms}}),\\ &p(\eta_0) = \text {N}\left(\mu_{\eta_0},\sigma^2_{\eta_0}\right),\\ &p(\lambda_{0,main})~=~ \text {N}\left(\mu_{\lambda_{0,main}},\sigma^2_{\lambda_{0,main}}\right)\mathcal{I}\left(c,\infty\right),\\ &p(\lambda_{0,inter})~=~ \text {N}\left(\mu_{\lambda_{0,inter}},\sigma^2_{\lambda_{0,inter}}\right),\\ \end{aligned} \end{align} $$

where $\textbf {I}_{D+1}$ is a $(D+1)$ -dimensional identity matrix. Parameter c is a truncation parameter. Some literature restricts the main effect terms of $\boldsymbol {\lambda }$ to non-negative values (Zhan et al., Reference Zhan, Jiao, Man and Wang2019). To address this, a truncation parameter c is introduced to adjust the range of values for the prior parameter $\lambda _{0,main}$ . For example, when c is set to $-\infty $ , there is no restriction on $\lambda _{0,main}$ , while setting $c=0$ restricts $\lambda _{0,main}$ to non-negative values. In practice, users can adjust the value of c to restrict the range of $\lambda _{0,main}$ according to their specific requirements. Let $\boldsymbol {\Omega }=(\boldsymbol {\delta }_0,\mu _{\eta },\sigma ^2_{\eta },\mu _{\lambda },\sigma ^2_{\lambda })$ , the joint posterior distribution of $(\textbf {Y},\boldsymbol {z},\boldsymbol {\pi },\boldsymbol {\lambda }^*,\boldsymbol {\lambda }^*_{0} | \textbf {Q},\boldsymbol {\Omega },\widetilde {\textbf {A}})$ based on the tight lower bound can be represented by

(20) $$ \begin{align} & \underline{p(\textbf{Y},\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}^*_{0} | \textbf{Q},\boldsymbol{\Omega},\widetilde{\textbf{A}})}~ = ~ \underline{p(\textbf{Y}|\boldsymbol{z},\widetilde{\textbf{A}},\boldsymbol{\lambda}^*,\textbf{Q})} p(\boldsymbol{z}|\boldsymbol{\pi}) p(\boldsymbol{\pi}) p(\boldsymbol{\lambda}^*|\boldsymbol{\lambda}^*_{0})p(\boldsymbol{\lambda}^*_{0})\nonumber\\ &~ \propto~\prod_{i=1}^N\prod_{j=1}^J \prod_{l=1}^L \left\{\sigma(\xi_{ijl})\exp\left(\frac{y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*-\xi_{ijl}}{2}-\tau(\xi_{ijl})(\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*\boldsymbol{h}_{jl}^{*\text{T}}\boldsymbol{\lambda}_j^*-\xi_{ijl}^{2})\right)\right\}^{z_{il}}\nonumber\\ &\times \prod_{i=1}^N\prod_{l=1}^L\pi_l^{z_{il}}\prod_{l=1}^L\pi_l^{\delta_{0l}}\prod_{j=1}^J\exp\left\{-\frac{(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)^{\text{T}}(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)}{2}\right\} \exp\left\{-\frac{(\eta_{0}-\mu_{\eta_0})^2}{2\sigma^2_{\eta_0}}\right\}\nonumber\\ &\times \exp\left\{-\frac{(\lambda_{0,main}-\mu_{\lambda_{0,main}})^2}{2\sigma^2_{\lambda_{0,main}}}\right\}\mathcal{I}\left(\lambda_{0,main}> c\right) \exp\left\{-\frac{(\lambda_{0,inter}-\mu_{\lambda_{0,inter}})^2}{2\sigma^2_{\lambda_{0,inter}}}\right\}\times const, \end{align} $$

where $const$ denotes a constant. The logarithm of $\underline {p(\textbf {Y},\boldsymbol {z},\boldsymbol {\pi },\boldsymbol {\lambda }^*,\boldsymbol {\lambda }^*_{0} | \textbf {Q},\boldsymbol {\Omega },\widetilde {\textbf {A}})} $ can be further expressed as

(21) $$ \begin{align} \begin{aligned} & \log \underline{p(\textbf{Y},\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}^*_{0} | \textbf{Q},\boldsymbol{\Omega},\widetilde{\textbf{A}})}\\ &= \sum_{i=1}^N\sum_{j=1}^J\sum_{l=1}^L{z_{il}}\left\{ \log(\sigma(\xi_{ijl}))+\frac{y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*-\xi_{ijl}}{2}-\tau(\xi_{ijl})(\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*\boldsymbol{h}_{jl}^{*\text{T}}\boldsymbol{\lambda}_j^*-\xi_{ijl}^{2})\right\}\\ &\quad+ \sum_{i=1}^N\sum_{l=1}^Lz_{il}\log \pi_l+\sum_{l=1}^L\delta_{0l}\log \pi_l -\sum_{j=1}^J\frac{(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)^{\text{T}}(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)}{2}\\ &\quad-\frac{(\eta_{0}-\mu_{\eta_0})^2}{2\sigma^2_{\eta_0}}- \frac{(\lambda_{0,main}-\mu_{\lambda_{0,main}})^2}{2\sigma^2_{\lambda_{0,main}}}\mathcal{I}\left(\lambda_{0,main}>c\right) -\frac{(\lambda_{0,inter}-\mu_{\lambda_{0,inter}})^2}{2\sigma^2_{\lambda_{0,inter}}}+const.\\ \end{aligned} \end{align} $$

3.5 Implementation of VBEM-M algorithm for LCDM

Assuming that the joint variational density of $(\boldsymbol {z},\boldsymbol {\pi },\boldsymbol {\lambda }^*,\boldsymbol {\lambda }^*_0)$ for the LCDM satisfies mean-field theory, the following equation holds:

(22) $$ \begin{align} q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}^*_0)=\left(\prod_{i=1}^N q(\boldsymbol{z})\right) q(\boldsymbol{\pi}) \left (\prod_{j=1}^{J}q(\boldsymbol{\lambda}^*_j)\right)q(\eta_0)q(\lambda_{0,main}) q(\lambda_{0,inter}). \end{align} $$

Let $\varTheta =(\boldsymbol {z},\boldsymbol {\pi },\boldsymbol {\lambda }^*,\boldsymbol {\lambda }^*_0)$ ; in terms of Eqs. (9) and (21), the ELBO $\underline {\mathcal {L}}(q(\varTheta ))$ can then be derived as

(23) $$ \begin{align} \begin{aligned} &\underline{\mathcal{L}}(q(\Theta))=\text{E}_{q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}_0^*)}[\log p(\textbf{Y},\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}^*_{0} | \textbf{Q},\boldsymbol{\Omega},\widetilde{\textbf{A}})-\log q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}_0^*)]\\ &\geq \text{E}_{q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}_0^*)}[\log \underline{p(\textbf{Y},\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}^*_{0} | \textbf{Q},\boldsymbol{\Omega},\widetilde{\textbf{A}})}-\log q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}_0^*)]\\ &= \text{E}_{q(\boldsymbol{z},\boldsymbol{\pi},\boldsymbol{\lambda}^*,\boldsymbol{\lambda}_0^*)}\Big[\sum_{i=1}^N\sum_{j=1}^J\sum_{l=1}^L z_{il}\left\{ \log(\sigma(\xi_{ijl}))+\frac{y_{ij}\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*-\xi_{ijl}}{2}-\tau(\xi_{ijl})(\boldsymbol{\lambda}_j^{*\text{T}}\boldsymbol{h}_{jl}^*\boldsymbol{h}_{jl}^{*\text{T}}\boldsymbol{\lambda}_j^*-\xi_{ijl}^{2})\right\} \\ &\quad+\sum_{i=1}^N\sum_{l=1}^Lz_{il}\log \pi_l+\sum_{l=1}^L(\delta_{0l}-1)\log \pi_l -\sum_{j=1}^J\frac{(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)^{\text{T}}(\boldsymbol{\lambda}^*_j-\boldsymbol{\lambda}_{0}^*)}{2}\\ &\quad-\frac{(\eta_{0}-\mu_{\eta_0})^2}{2\sigma^2_{\eta_0}}-\frac{(\lambda_{0,main}-\mu_{\lambda_{0,main}})^2}{2\sigma^2_{\lambda_{0,main}}}\mathcal{I}\left(\lambda_{0,main}> c\right)-\frac{(\lambda_{0,inter}-\mu_{\lambda_{0,inter}})^2}{2\sigma^2_{\lambda_{0,inter}}}\\ &\quad-\sum_{i=1}^N\sum_{l=1}^Lz_{il}\log \pi^*_{il} -\sum_{l=1}^L(\delta^*_{l}-1)\log \pi_l -\sum_{j=1}^J\frac{(\boldsymbol{\lambda}^*_j-\boldsymbol{m}_j^*)^{\text{T}}\boldsymbol{V}_j^{*-1}(\boldsymbol{\lambda}^*_j-\boldsymbol{m}_{j}^*)}{2}\\ &\quad+\frac{(\eta_{0}-\mu^*_{\eta_0})^2}{2\sigma^{*2}_{\eta_0}}+\frac{(\lambda_{0,main}-\mu^*_{\lambda_{0,main}})^2}{2\sigma^{*2}_{\lambda_{0,main}}}\mathcal{I}\left(\lambda_{0,main}> c\right)-\frac{(\lambda_{0,inter}-\mu^*_{\lambda_{0,inter}})^2}{2\sigma^{*2}_{\lambda_{0,inter}}}\Big]+const\\ &\triangleq \underline{\mathcal{L}}^*(q(\Theta),\boldsymbol{\xi}), \end{aligned} \end{align} $$

where $\underline {\mathcal {L}}^*(q(\Theta ),\boldsymbol {\xi })$ is a tight lower bound of $\underline {\mathcal {L}}(q(\Theta ))$ . Next, we maximize $\underline {\mathcal {L}}^*(q(\Theta ),\boldsymbol {\xi })$ to obtain estimates of latent variables $\boldsymbol {z}$ , model parameters $(\boldsymbol {\pi }, \boldsymbol {\lambda }^*)$ , hyperparameters $\boldsymbol {\lambda }_0^*$ , and local point parameter $\boldsymbol {\xi }$ . Specifically, there are three steps to the implementation process:

1. (a) VBE step: update variational density for latent variable;

2. (b) VBM step: update variational densities for model parameters and hyperparameters;

3. (c) M step: update local point parameter $\boldsymbol {\xi }$ by maximizing $\underline {\mathcal {L}}^*(q(\Theta ),\boldsymbol {\xi })$ .

In the following text, $q^*(\cdot )$ denotes the optimal variational posterior in each iteration. To keep things simple, we only present the core formulation for updating. The specifics can be found in the Supplementary Material. The estimation procedure of the VBEM-M algorithm is shown in Table 1. In Table 1 and subsequent tables, all parameters are estimated using their posterior means. In addition, the specific implementation process of each step for the VBEM-M algorithm is shown in Figure 1.

Table 1 Estimation procedure of the VBEM-M algorithm

Figure 1 Graphical illustration of the VBEM-M algorithm implementation process. Let $\Theta ^{*(t)}=(\boldsymbol {\pi }^{(t)}, \boldsymbol {\lambda }^{*(t)}, \boldsymbol {\lambda }_0^{*(t)})$ . the variational density of the latent variable $\boldsymbol {z}^{(t+1)}$ is updated in the VBE-step. In VBM-step, the variational densities for model parameters and hyperparameters $\Theta ^{*(t+1)}$ are updated. In M-step, we update $\boldsymbol {\xi }^{*(t+1)}$ by maximizing $\text {L}(q(\boldsymbol {z}^{(t+1)},\Theta ^{*(t+1)}),\boldsymbol {\xi }^{*(t)})$ .

(a) VBE step. In this step, we update the variational density of $\boldsymbol {z}_i$ for each i, where $i=1,\cdots ,N$ . $q^*(\boldsymbol {z}_i)$ is derived to be a categorical distribution with parameter $\boldsymbol {\pi }^*_i$ . That is,

(24) $$ \begin{align} q^*(\boldsymbol{z}_i|\boldsymbol{\pi}^*_i)=\prod_{l=1}^L\pi_{il}^{*z_{il}}, \end{align} $$

where

(25) $$ \begin{align} \begin{aligned} &\qquad \qquad \quad \pi_{il}^*=\frac{\rho_{il}}{\sum_{l=1}^{L}\rho_{il}},\\ \rho_{il} =\exp\Bigg\{ & \sum_{j=1}^J\Big\{ \log(\sigma(\xi_{ijl}))+\frac{y_{ij}\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*]^{\text{T}}\boldsymbol{h}_{jl}^*-\xi_{ijl}}{2}\\ &-\tau(\xi_{ijl})(\boldsymbol{h}_{jl}^{*\text{T}}\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*\boldsymbol{\lambda}_j^{*\text{T}}]\boldsymbol{h}_{jl}^*-\xi_{ijl}^{2})\Big\}+\text{E}_{q(\boldsymbol{\pi})}\log (\pi_l)\Bigg\}. \end{aligned} \end{align} $$

(b) VBM step. In this step, we update the variational density for $\boldsymbol {\pi }$ , $\boldsymbol {\lambda }^*_j\,(j=1,\cdots ,J)$ , $\eta _0$ , $\lambda _{0,main}$ and $\lambda _{0,inter}.$

(b1) Update the variational density for  $\boldsymbol {\pi }$

$q^*(\boldsymbol {\pi })$ is derived to be a Dirichlet distribution with parameter $\boldsymbol {\delta }^*$ . That is,

(26) $$ \begin{align} q^*(\boldsymbol{\pi}|\boldsymbol{\delta}^*)\propto\prod_{l=1}^L\pi_l^{\delta_{l}^*-1}, \end{align} $$

where

(27) $$ \begin{align} \delta_l^* = \sum_{i=1}^{N}\text{E}_{q(\boldsymbol{z}_{i})}[z_{il}]+\delta_{0l}. \end{align} $$

(b2) Update the variational density for  $\boldsymbol {\lambda }^*_j$

$q(\boldsymbol {\lambda }^*_j)$ is proportional to a multivariate normal distribution with mean vector $\boldsymbol {m}_{j}^*$ and covariance $\boldsymbol {V}_j^*$ . That is,

(28) $$ \begin{align} q^*(\boldsymbol{\lambda}^*_j|\boldsymbol{m}_{j}^*,\boldsymbol{V}_j^*)\propto\exp\left\{-\frac{(\boldsymbol{\lambda}^*_j-\boldsymbol{m}_{j}^*)^{\text{T}}\boldsymbol{V}_j^{*-1}(\boldsymbol{\lambda}^*_j-\boldsymbol{m}_{j}^*)}{2}\right\}, \end{align} $$

where

(29) $$ \begin{align} \begin{aligned} & \boldsymbol{V}_j^{*-1} =\textbf{I}_{D+1}^{-1}+2 \sum_{i=1}^{N}\sum_{l=1}^{L} \text{E}_{q(\boldsymbol{z}_{i})}[z_{il}]\tau(\xi_{ijl})\boldsymbol{h}_{jl}^*\boldsymbol{h}_{jl}^{*\text{T}},\\ &\boldsymbol{m}^*_j= \boldsymbol{V}_j^{*}\left (\boldsymbol{\lambda}^*_{0}+\frac{1}{2}\sum_{i=1}^{N}\sum_{l=1}^{L}\text{E}_{q(\boldsymbol{z}_{i})}[z_{il}]y_{ij}\boldsymbol{h}_{jl}^*\right ). \end{aligned} \end{align} $$

(b3) Update the variational density for  $\mathbf {\eta _0}$

$q^*(\eta _0)$ is proportional to a normal distribution with mean $\mu _{\eta _0}^*$ and variance $\sigma _{\eta _0}^{*2}$ . That is,

(30) $$ \begin{align} q^*(\eta_0|\mu_{\eta_0}^*,\sigma_{\eta_0}^{*2})\propto\exp\left\{-\frac{(\eta_0-\mu_{\eta_0}^*)^2}{2\sigma_{\eta_0}^{*2}}\right\}, \end{align} $$

where

(31) $$ \begin{align} \begin{aligned} &(\sigma^{*2}_{\eta_0})^{-1}=J+\frac{1}{\sigma^2_{\eta_0}},\\ &\mu^*_{\eta_0}= \sigma^{*2}_{\eta_0}\left(\frac{\mu_{\eta_0}}{{\sigma}^2_{\eta_0}}+\sum_{j=1}^{J}\left(\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*]\right)_{\eta}\right), \end{aligned} \end{align} $$

where $\text {E}_{q(\boldsymbol {\lambda }_j^*)}[\boldsymbol {\lambda }_j^*]_{\eta }$ is the corresponding expected value of the element $\eta _j$ in the vector $\boldsymbol {\lambda }^*$ .

(b4) Update the variational density for  $\mathbf {\lambda _{0,main}}$

$q^*(\lambda _{0,main})$ is proportional to a truncated normal distribution with mean $\mu _{\lambda _{0,main}}^*$ and variance $\sigma _{\lambda _{0,main}}^{*2}$ . Specifically,

(32) $$ \begin{align} q^*(\lambda_{0,main}|\mu_{\lambda_{0,main}}^*,\sigma_{\lambda_{0,main}}^{*2})\propto\exp\left\{-\frac{(\lambda_{0,main}-\mu_{\lambda_{0,main}}^*)^2}{2\sigma_{\lambda_{0,main}}^{*2}}\right\}\mathcal{I}\left(c,\infty\right), \end{align} $$

where

(33) $$ \begin{align} \begin{aligned} &(\sigma^{*2}_{\lambda_{0,main}})^{-1}=J^*_{main}+\frac{1}{\sigma^2_{\lambda_{0,main}}},\\ &\mu^*_{\lambda_{0,main}}= \sigma^{*2}_{\lambda_{0,main}}\left(\frac{\mu_{\lambda_{0,main}}}{{\sigma}^2_{\lambda_{0,main}}}+\sum_{d=1}^{K}\sum_{j\in J_d}\left(\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*]\right)_{\lambda_d}\right), \end{aligned} \end{align} $$

where $J^*_{main}$ denotes the number of all main effect terms, $J_d=\{j: \lambda _{jd}\neq 0 \}$ , and $\text {E}_{q(\boldsymbol {\lambda }_j^*)}[\boldsymbol {\lambda }_j^*]_{\lambda _d}$ is the corresponding expected value of the element $\lambda _{jd}$ in the vector $\boldsymbol {\lambda }^*$ .

(b5) Update the variational density for  $\mathbf {\lambda _{0,inter}}$

$q^*(\lambda _{0,inter})$ is proportional to a truncated normal distribution with mean $\mu _{\lambda _{0,inter}}^*$ and variance $\sigma _{\lambda _{0,inter}}^{*2}$ . Specifically,

(34) $$ \begin{align} q^*(\lambda_{0,inter}|\mu_{\lambda_{0,inter}}^*,\sigma_{\lambda_{0,inter}}^{*2})\propto\exp\left\{-\frac{(\lambda_{0,inter}-\mu_{\lambda_{0,inter}}^*)^2}{2\sigma_{\lambda_{0,inter}}^{*2}}\right\}, \end{align} $$

where

(35) $$ \begin{align} \begin{aligned} &(\sigma^{*2}_{\lambda_{0,inter}})^{-1}=J^*_{inter}+\frac{1}{\sigma^2_{\lambda_{0,inter}}},\\ &\mu^*_{\lambda_{0,inter}}= \sigma^{*2}_{\lambda_{0,inter}}\left(\frac{\mu_{\lambda_{0,inter}}}{{\sigma}^2_{\lambda_{0,inter}}}+\sum_{d=K+1}^{D}\sum_{j\in J_d}\left(\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*]\right)_{\lambda_d}\right), \end{aligned} \end{align} $$

where $J^*_{inter}$ denotes the number of all interaction terms.

(c) M step. In this step, we update the local point parameter $\xi _{ijl}\,~(i=1,\cdots ,N;\;j=1,\cdots ,J;\;l=1,\cdots ,L)$ . To obtain the optimal $\xi _{ijl}$ , we need to maximize $\underline {\mathcal {L}}^*(q(\Theta ),\boldsymbol {\xi })$ by computing the derivative of $\xi _{ijl}$ to zero:

(36) $$ \begin{align} \frac{\partial\underline{\mathcal{L}}^*(q(\Theta),\boldsymbol{\xi})}{\partial\xi_{ijl}}=0. \end{align} $$

Therefore, we have

(37) $$ \begin{align} \xi_{ijl}^2=\xi_{jl}^2=\boldsymbol{h}_{jl}^{*\text{T}} \text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*\boldsymbol{\lambda}_j^{*\text{T}}]\boldsymbol{h}_{jl}^*. \end{align} $$

Considering the aforementioned presentation of the VBEM-M algorithm, it is clear that we need to compute a large number of expectations using categorical, Dirichlet, normal, multivariate normal, and truncated normal distributions. Some formulae for calculating these expectations are as follows:

(38) $$ \begin{align} \begin{aligned} &\text{E}_{q(\boldsymbol{z}_{i})}[z_{il}]=\pi^*_{il}, ~~\text{E}_{q(\boldsymbol{\pi})}\log (\pi_l)=\psi(\delta_l^*)-\psi\left(\sum_{l=1}^{L}\delta_l^*\right),\\ &\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*]=\boldsymbol{\mu}_{j}^*,~~\text{E}_{q(\boldsymbol{\lambda}_j^*)}[\boldsymbol{\lambda}_j^*\boldsymbol{\lambda}_j^{*\text{T}}]=\boldsymbol{\Sigma}_j^*+\boldsymbol{\mu}_j^*\boldsymbol{\mu}_j^{*\text{T}},\\ &\text{E}_{q(\eta_0)}[\eta_0]=\mu_{\eta_0}^*,~~\text{E}_{q(\eta_0)}[(\eta_0-\mu_{\eta_0}^*)^2]=\sigma^{*2}_{\eta_0},\\ &\text{E}_{q(\lambda_{0,main})}[\lambda_{0,main}]=\mu_{\lambda_{0,main}}^*+\sigma^{*2}_{\lambda_{0,main}}\,\frac{\phi(u)}{\Phi(u)},u=\frac{c-\mu_{\lambda_{0,main}}^*}{\sigma^{*}_{\lambda_{0,main}}},\\ &\text{E}_{q(\lambda_{0,main})}[(\lambda_{0,main}-\mu_{\lambda_{0,main}}^*)^2]=\sigma^{*2}_{\lambda_{0,main}}\left(1-\mu_{\lambda_{0,main}}^*\frac{\phi(u)}{\Phi(u)}\right),\\ &\text{E}_{q(\lambda_{0,inter})}[\lambda_{0,inter}]=\mu_{\lambda_{0,inter}}^*,~~\text{E}_{q(\lambda_{0,inter})}[(\lambda_{0,inter}-\mu_{\lambda_{0,inter}}^*)^2]=\sigma^{*2}_{\lambda_{0,inter}}, \end{aligned} \end{align} $$

where $\psi (\cdot )$ is $\psi (x)=\frac {\mathrm {d}}{\mathrm {d}x}\log \Gamma (x)$ , $\Gamma (x)=\int _0^{\infty } t^{(x-1)}\exp (-t)\, dt $ , $\phi (\cdot )$ is the density function of a standard normal distribution, and $\Phi (\cdot )$ is the cumulative distribution function of a standard normal distribution.

4.1 Data generation

Item response data $x_{ij}$ is generated from a Bernoulli distribution with probability of correct response $P(x_{ij}=1|\boldsymbol {\alpha }_i,\boldsymbol {\lambda }^{*},\boldsymbol {q}_j)$ . The true values of the item parameters based on DINA model are constrained by considering four different levels of noise to investigate the correlation between noise and recovery. For each item, the following scenarios are considered. (a1) Low noise level (LNL): $s_j=g_j=0.1$ , with corresponding true values $\eta _j= -2.1972$ , $\lambda _j=4.3944$ . (a2) High noise level (HNL): $s_j=g_j=0.2$ , with corresponding true values $\eta _j= -1.3863$ , $\lambda _j=2.7726$ . (a3) Slipping higher than guessing (SHG): $s_j=0.2$ , $g_j=0.1$ , with corresponding true values $\eta _j= -2.1972$ , $\lambda _j=3.5835$ . (a4) Guessing higher than slipping (GHS): $s_j=0.1$ , $g_j=0.2$ , with corresponding true values $\eta _j= -1.3863$ , $\lambda _j=3.5835$ .

To generate the attribute-mastery patterns, we used the same procedure as Chiu and Douglas (Reference Chiu and Douglas2013), which takes into account the correlations among the attributes. Specifically, $\boldsymbol {\alpha }^*_{i}=(\alpha ^*_{i1},\cdots ,\alpha ^*_{ik},\cdots ,\alpha ^*_{iK})^{\text {T}}$ are generated from a multivariate normal distribution; that is, $\boldsymbol {\alpha }^*_i\sim \text {N}(\boldsymbol {0}_K,\boldsymbol {\Sigma }_{K\times K})$ , where $\boldsymbol {0_K}=(0,\dots ,0)^{\text {T}}_{K\times 1}$ and

$$ \begin{align*}\boldsymbol{\Sigma_{K\times K}}=\begin{bmatrix} 1 & \dots & \sigma\\ \vdots& \ddots & \vdots\\ \sigma & \dots& 1 \end{bmatrix}_{K \times K}, \end{align*} $$

where the off-diagonal elements of $\Sigma _{K\times K}$ are $\sigma $ . As $\sigma $ increases from 0 to 1, the correlation between attributes also increases from 0 to maximum. The relationships between the attribute profiles $\boldsymbol {\alpha }_i$ and $\boldsymbol {\alpha }^*_i$ can be expressed as $\alpha _{ik}=1$ if $\alpha _{ik}^*>0$ and $\alpha _{ik}=0$ otherwise. Although the $\text {Q}$ -matrices are created randomly, they still conform to the identifiability constraints outlined by Chen et al. (Reference Chen, Liu, Xu and Ying2015, Reference Chen, Culpepper, Chen and Douglas2017), Liu and Andersson (Reference Liu, Andersson and Skrondal2020), and Xu and Shang (Reference Xu and Shang2018). We present the $\text {Q}$ -matrices used in these simulations in the Supplementary Material.

4.2 Prior distributions

The prior parameter $\boldsymbol {\delta }_0$ is set as $\boldsymbol {\delta }_0=\boldsymbol {1}_L$ (Culpepper, Reference Culpepper2015; Zhan et al., Reference Zhan, Jiao, Man and Wang2019), where $\boldsymbol {1}_L$ denotes a L-dimensional vector with all elements equal to 1. The hyperparameters are chosen as follows: $\mu _{\eta _0}=-2$ , $\mu _{\lambda _{main}}=\mu _{\lambda _{inter}}=0$ , and $\sigma ^2_{\eta _0}=\sigma ^2_{\lambda _{main}}=\sigma ^2_{\lambda _{inter}}=10$ .

4.3 Estimation software

We implemented four different approaches, namely, the VBEM-M algorithm, VB algorithm, MCMC sampling algorithm, and EM algorithm, using the R programming language (R Core Team, 2017) on a desktop computer equipped with Intel (R) Core (TM) i5-10400 CPU @ 2.90GHz, 16GB RAM. To enhance the computational efficiency of the VBEM-M method, we utilized two R packages, “Rcpp” (Eddelbuettel & Francois, Reference Eddelbuettel and Francois2011) and “RcppArmadillo” (Eddelbuettel & Sanderson, Reference Eddelbuettel and Sanderson2014), to call the C++ programming language. The R code of our VBEM-M algorithm can be found in the Supplementary Material. We used the R package “variationalDCM” (Hijikata et al., Reference Hijikata, Oka, Yamaguchi and Okada2023) to implement the VB method. The MCMC sampling algorithms were implemented separately using the R packages “dina” (Culpepper & Balamuta, Reference Culpepper and Balamuta2019), which is integrated with the C++ program, and “R2jags” (Su & Yajima, Reference Su and Yajima2015) which is associated with the JAGS program (Plummer, Reference Plummer2003). The EM algorithm was implemented using the R packages “GDINA” (Ma & de la Torre, Reference Ma and de la Torre2020) and “CDM” (George et al., Reference George, Robitzsch, Kiefer, Groß and Ünlü2016), respectively.

4.4 Convergence diagnosis

The VBEM-M algorithm was considered converged if the absolute difference between two consecutive iterations was less than $e_0=10^{-4}$ , or if the number of iterations T had reached 2,000. When using the R packages “dina” and “R2jags” to implement the MCMC sampling algorithms, for the DINA model, Culpepper (Reference Culpepper2015) demonstrated that it would have converged after 750 iterations, thus the chain length was set to 2,000 and the first 1,000 iterations were set as a “burn-in” period. For the saturated LCDM, we chose a chain length of 10,000, with a burn-in of 5,000. For the EM algorithm, when employing the R package “GDINA,” the convergence criteria is when the maximum absolute change in item success probabilities between consecutive iterations was smaller than $e_0=10^{-4}$ or when T exceeded 2,000. In addition, when using the R package “CDM,” iteration will end if the maximal change in parameter estimates is below $e_0=0.001$ .

4.5 Evaluation Criteria

For item parameters and class membership probability parameters, we assess the accuracy of parameter estimation using bias and root mean square error (RMSE). For attribute parameters, we adopt the following two evaluation indices: the pattern-wise agreement rate (PAR), which indicates the rates of correct classification for attribute patterns, and is formulated as

(39) $$ \begin{align} \text{PAR}=\frac{1}{N}\sum_{i=1}^{N}\mathcal{I}(\hat{\boldsymbol{\alpha}}_i=\boldsymbol{\alpha}_i), \end{align} $$

and the attribute-wise agreement rate (AAR), which signifies the rates of correct classification for individual attributes, and is defined as

(40) $$ \begin{align} \text{AAR}(k)=\frac{1}{N}\sum_{i=1}^{N}\mathcal{I}(\hat{\alpha}_{ik}=\alpha_{ik}), \end{align} $$

where $\boldsymbol {\alpha }_i$ is the true value of the ith student’s attribute profile and $\hat {\boldsymbol {\alpha }}_i$ is the estimated value of $\boldsymbol {\alpha }_i$ . $\hat {\alpha }_{ik}$ is the estimated value of $\alpha _{ik}$ for the specific attribute k.

4.6 Simulation study 1

In this simulation study, we explored the performance of the VBEM-M algorithm under various simulation conditions. We set the test length to $J=30$ , the number of attributes was set to $K=5$ , and the corresponding $\text {Q}$ -matrix is shown in the Supplementary Material. The following manipulated conditions were considered: (A) number of examinees $N=1,000$ and 2,000; (B) correlation among attributes $\sigma =0$ , 0.3 and 0.7; and (C) noise levels LNL, HNL, SHG, and GHS. Fully crossing different levels of these three factors yields 24 conditions (2 sample sizes $\times $ 3 correlations $\times $ 4 noise levels). There were 100 replications for each simulation condition. The recovery results of parameters are displayed in Tables 2 and 3 and Figure 3.

Table 2 The accuracy of item parameters and class membership probability parameters using the VBEM-M algorithm in simulation study 1

Note: The values outside parentheses represent the RMSE, while the values inside the parentheses indicate the bias. These reflect the average RMSE, bias and SD for all intercept parameters $\boldsymbol {\eta }$ , slope parameters $\boldsymbol {\lambda }$ , and class membership probability parameters $\boldsymbol {\pi }$ .

Table 3 The accuracy of attribute profile parameters using the VBEM-M algorithm in simulation study 1

Note: AAR1 represents the correct classification rate for the first attribute, AAR2 for the second attribute, AAR3 for the third attribute, AAR4 for the fourth attribute, and AAR5 for the fifth attribute. PAR stands for the pattern-wise agreement rate.

Figure 2 The bias and RMSE of $\eta $ and $\lambda $ for each item in the simulation study 1. The $\text {Q}$ -Matrix denotes the skills required for each item along the x axis, where the black square =“1” and white square =“0”.

Figure 3 The boxplots of bias and RMSE for $\eta $ , $\lambda $ and $\pi $ estimated by the VBEM-M, VB, MCMC-dina, MCMC-R2jags, EM-GDINA and EM-CDM with $\sigma =0.3$ under the LNL condition in simulaion study 2.

The following conclusions can be drawn from Tables 2 and 3. (1) Given the correlation and noise levels, when the number of examinees is increased from 1,000 to 2,000, the average RMSE, the average bias, and standard deviation (SD) for η , λ , and π show decreasing trends. For example, when the correlation among attributes is 0.3 and the LNL is applied, increasing the number of examinees from 1,000 to 2,000 results in the average bias of η decreasing from -0.0140 to -0.0077, and the average bias of λ decreasing from 0.0307 to 0.0133. The average RMSE of η decreases from 0.1369 to 0.0981, the average RMSE of λ from 0.2337 to 0.1669, and the average RMSE of π from 0.0022 to 0.0016. The SD of η decreases from 0.0937 to 0.0664, the SD of λ decreases from 0.1617 to 0.1152, and the SD of π decreases from 0.0051 to 0.0037. (2) When the number of examinees and the noise level are given, with increasing $\sigma $ , the average RMSE for $\boldsymbol {\eta }$ increase somewhat. This indicates that $\boldsymbol {\eta }$ is less impacted by the correlation between attributes. $\boldsymbol {\lambda }$ is substantially more impacted by $\sigma $ ; specifically, the average bias and RMSE for $\boldsymbol {\lambda }$ tend to decrease markedly as $\boldsymbol {\sigma }$ increases. In the meanwhile, RMSE for $\boldsymbol {\pi }$ also tend to decrease as $\boldsymbol {\sigma }$ increases. For example, when the number of examinees is fixed at 1,000 and the LNL noise level is applied, the average bias are –0.0118, 0.140, 0.104, respectively, and the average RMSE rises from 0.1351 to 0.1388 when $\sigma $ increases. The change in bias and RMSE of $\boldsymbol {\eta }$ are found to be slight. However, the decreases in bias and RMSE are markedly greater for $\boldsymbol {\lambda }$ , with the average bias of $\boldsymbol {\lambda }$ decreasing from 0.0365 to 0.0279 and the corresponding average RMSE decreasing from 0.2560 to 0.2216. For $\boldsymbol {\pi }$ , the average bias remains at 0.0000 in all conditions, while the average RMSE exhibits the largest change in the HNL condition, decreasing from 0.0058 to 0.0048. (3) The accuracy of attribute profile recovery is highest under the LNL condition because the noise is the lowest. For example, with a fixed number of examinees at 1,000 and a correlation of $\sigma =0$ , the PAR is 0.9025 under the LNL condition and only 0.6736 under the HNL condition. Under the LNL condition, the AAR values for five attributes exceed 0.9667 across various sample sizes and levels of attribute correlation. Moreover, the accuracy of attribute profile recovery tends to improve as $\sigma $ increases.

In Figure 3, as an explanation, we only show the recovery results for the LNL and HNL based on the sample size $N=1,000$ . On each item, the bias of $\boldsymbol {\eta }$ are almost the same for the LNL and the HNL. Furthermore, when the correlation between attributes is strengthened ( $\sigma $ from 0 to 0.7), there is no difference between the bias and RMSE of $\boldsymbol {\eta }$ in the LNL (HNL). It was also discovered that, for both low and high levels of noise, the RMSE of $\boldsymbol {\eta }$ is lower when the items evaluate more attributes. At low noise levels, for instance, the RMSE of $\boldsymbol {\eta }$ for the first item evaluating one attribute is greater than that for the eleventh item evaluating the first three attributes together. For $\boldsymbol {\lambda }$ , although the bias of $\boldsymbol {\lambda }$ differs on each item at low and high noise levels, the values of bias are basically around 0. Similarly, for both low and high levels of noise, the RMSE of $\boldsymbol {\lambda }$ is lower when items have higher correlation among themselves. This is because as the attribute correlation increases, more accurate estimates of ${\alpha }$ are obtained, which in turn enhances the accuracy of ${\lambda }$ estimates. This also provides an empirical guarantee for our later practical research. That is, when designing the items, we should aim to achieve higher correlations between attributes to increase the accuracy of parameter estimation.

Additionally, we assess the performance of the VBEM-M algorithm under different initial values (please see the Supplementary Material for details), and the results showed that our VBEM-M algorithm is not affected by the different initial values.

4.7 Simulation study 2

The purpose of this simulation study is to compare the proposed method with Yamaguchi and Okada’s (Reference Yamaguchi and Okada2020b) VB method, the MCMC sampling algorithms, and the EM algorithm in terms of parameter accuracy for the DINA model. Specifically, the R package “variationalDCM” was used to implement Yamaguchi and Okada’s (Reference Yamaguchi and Okada2020b) VB method, while the R packages “dina” and “R2jags” were used to implement the MCMC sampling algorithms. The EM algorithm was implemented using the R packages “GDINA” and “CDM”.

The simulation design is as follows: the test length was fixed at $J=30$ , and the number of attributes was set to $K=5$ . The varying conditions of the simulation are as follows: (D) The number of examinees $N=200$ , 500, 1,000, and 2,000; (E) correlation among attributes $\sigma =0$ , 0.3, and 0.7; and (F) LNL and HNL conditions. Fully crossing different levels of these two factors yields 24 conditions (4 sample sizes $\times $ 3 correlations $\times $ 2 noise levels). Each simulation condition was replicated 100 times. The recovery results of item parameters and attribute profile recovery for all six methods are shown in Tables 4 and 5. Due to the space limit, we only present the results with the correlation $\sigma =0.3$ in Tables 4 and 5; the other two correlation cases ( $\sigma =0$ and $\sigma =0.7$ ) are given in the Supplementary Material. Figure 2 depicts the boxplots of the bias and RMSE for $\boldsymbol {\eta }$ , $\boldsymbol {\lambda }$ , and $\boldsymbol {\pi }$ estimated by the six methods with $\sigma =0.3$ under the LNL condition. Table 6 shows the computation time for these six methods under the same conditions. Here, the displayed computation time is the average time across 100 replications.

Table 4 The accuracy of item parameters and class membership probability parameters using the VBEM-M, VB, MCMC-dina, MCMC-R2jags, EM-GDINA, and EM-CDM algorithms for the DINA model under the $\sigma =0.3$ condition in simulation study 2

Note: The values outside the parentheses represent the RMSE, while the values inside the parentheses indicate the bias. Here, RMSE and Bias denote the average RMSE and Bias, respectively, for all intercept parameters $\boldsymbol {\eta }$ , all slope parameters $\boldsymbol {\lambda }$ and all class membership probability parameters $\boldsymbol {\pi }$ .

Table 5 The accuracy of attribute profile parameters using the VBEM-M, VB, MCMC-dina, MCMC-R2jags, EM-GDINA, and EM-CDM algorithms for the DINA model under the $\sigma =0.3$ condition in simulation study 2

Note: AAR1 represents the correct classification rate for the first attribute, AAR2 for the second attribute, AAR3 for the third attribute, AAR4 for the fourth attribute, and AAR5 for the fifth attribute. PAR stands for the pattern-wise agreement rate.

Table 6 The computational time (in seconds) for the VBEM-M, VB, MCMC-dina, MCMC-R2jags, EM-GDINA, and EM-CDM algorithms with the $\sigma =0.3$ condition based on DINA in simulation study 2

In Tables 4 and 5, as well as in the subsequent simulation studies, the RMSE and bias mentioned are the average RMSE and average bias. From Tables 4 and 5, we can draw the following conclusions: (1) The VBEM-M algorithm consistently outperforms the other five methods in terms of achieving lower RMSE values for item parameters $\boldsymbol {\eta }$ and $\boldsymbol {\lambda }$ under all four sample sizes, regardless of LNL or HNL condition. (2) For the EM algorithm, both EM-GDINA and EM-CDM methods have higher bias and RMSE for item parameters $\boldsymbol {\eta }$ and $\boldsymbol {\lambda }$ than four other methods, especially for a small sample size of N=200, under both LNL and HNL conditions. (3) With the same sample size and noise level, both MCMC methods (MCMC-dina and MCMC-R2jags) show similar estimation accuracy, as do the two EM methods (EM-GDINA and EM-CDM). (4) For parameter $\boldsymbol {\pi }$ , the estimated bias and RMSE of the six methods are basically the same under various identical simulation conditions, with no significant differences. (5) In terms of the accuracy of attribute profile recovery, the results of the six methods are essentially the same under each simulation condition.

From Table 6, we can see that the VBEM-M algorithm is highly efficient in terms of computation time. It performs faster than the VB method across most simulation conditions, and this speed advantage is more noticeable as sample sizes increase. Overall, the computation speed of the VBEM-M algorithm is second only to the two EM algorithms, i.e., EM-GDINA and EM-CDM. The two Bayesian methods, MCMC-dina and MCMC-R2jags, have longer computation time than the other four methods. Additionally, MCMC-dina is faster than MCMC-R2jags due to its use of the “Rcpp” and “RcppArmadillo” packages, which are built on C++ programming language.

4.8 Simulation study 3

This simulation study aims to evaluate the effectiveness of the VBEM-M algorithm on the saturated LCDM by comparing it with Yamaguchi and Okada’s (Reference Yamaguchi and Okada2020a) VB method, the MCMC sampling algorithms, and the EM algorithm. Specifically, the R package “variationalDCM” was used to implement Yamaguchi and Okada’s (Reference Yamaguchi and Okada2020a) VB method, the R package “R2jags” was used to implement the MCMC sampling algorithms, and the EM algorithm was implemented using the R package “GDINA”.

This simulation was designed with an attribute number of $K=3$ and a test length of $J=18$ . In the saturated LCDM, each item’s $\boldsymbol {\lambda }_j^*$ is an 8-dimensional vector ( $2^3=8$ ). The true values of $\boldsymbol {\lambda }^*$ are shown in Table 7. We conducted simulations across different sample sizes (N=1,000, 2,000) and attribute correlations ( $\sigma =0, 0.3, 0.7$ ), resulting in six different conditions. Each condition was replicated 100 times. Notably, an additional calculation procedure was needed for Yamaguchi and Okada’s (Reference Yamaguchi and Okada2020a) VB method, as the R package “variationalDCM” only reports the correct response probabilities for different attribute mastery patterns. We transformed these probabilities into LCDM parameters by solving a linear system of equations (Liu & Johnson, Reference Liu, Johnson, von Davier and Lee2019; Yamaguchi & Templin, Reference Yamaguchi and Templin2022). The parameter recovery results for the $\sigma =0.3$ condition are displayed in Tables 8 and 9. The estimation results for $\sigma =0$ and 0.7 are available in the Supplementary Material.

Table 7 True values of $\boldsymbol {\lambda }^*$ for the saturated LCDM in simulation study 3

Table 8 The accuracy of item parameters and class membership probability parameters using the VBEM-M, VB, MCMC-R2jags, and EM-GDINA algorithms for the LCDM model under the $\sigma =0.3$ condition in simulation study 3

Note: The values outside the parentheses represent the RMSE, while the values inside the parentheses indicate the bias. Here, RMSE and Bias denote the average RMSE and Bias, respectively, for all intercept parameters $\boldsymbol {\eta }$ , all main effect slope parameters $\boldsymbol {\lambda }_{main}$ , all interaction slope parameters $\boldsymbol {\lambda }_{inter}$ and all class membership probability parameters $\boldsymbol {\pi }$ .

Table 9 Evaluation of the accuracy of attribute profile parameters using the VBEM-M, VB, MCMC-R2jags and EM-GDINA Algorithms for the saturated LCDM under the $\sigma =0.3$ condition in simulation study 3

Note: AAR1 represents the correct classification rate for the first attribute, AAR2 for the second attribute, AAR3 for the third attribute. PAR stands for the pattern-wise agreement rate.

For a more detailed analysis, we split the parameter $\boldsymbol {\lambda }$ into two parts: $\boldsymbol {\lambda }_{main}$ (i.e., $\boldsymbol {\lambda }_1,\boldsymbol {\lambda }_2,\boldsymbol {\lambda }_3$ ) and $\boldsymbol {\lambda }_{inter}$ (i.e., $\boldsymbol {\lambda }_{12},\boldsymbol {\lambda }_{13},\boldsymbol {\lambda }_{23},\boldsymbol {\lambda }_{123}$ ), which represent the main effects and interactions, respectively. From the results, we can draw the following conclusions: (1) As the number of examinees increases, the RMSE for item parameters of all algorithms decreases. (2) The proposed VBEM-M algorithm performs better than other algorithms on all item parameters across all conditions, especially on the interactions. Specifically, in terms of the parameters $\boldsymbol {\eta }$ and $\boldsymbol {\pi }$ , VBEM-M has a slight advantage over the other algorithms, whereas it shows a significant advantage in estimating $\boldsymbol {\lambda }_{main}$ and $\boldsymbol {\lambda }_{inter}$ , particularly for the parameter $\boldsymbol {\lambda }_{inter}$ . On the other hand, the EM algorithm performs poorly with small sample sizes. For the $\boldsymbol {\lambda }_{inter}$ parameter, its RMSE exceeds 2 when $N=200$ . (3) Compared to other algorithms, VBEM-M performs significantly better with small sample sizes ${(}N=200, 500$ ), with noticeably lower RMSE. (4) It is worth noting that the results from all algorithms indicate that, although the interaction terms have smaller true values compared to the main effects, their estimation accuracy is worse. This suggests that estimating interaction effects is the most challenging aspect of the saturated LCDM model. (5) As for the accuracy of attribute profiles, there is no obvious difference among these algorithms, but VBEM-M still shows slightly higher accuracy than the others.

Table 10 shows the average computation time across 100 replications for the four algorithms under the $\sigma =0.3$ condition. The results indicate that our algorithm performs better than the other algorithms in terms of computational efficiency. Additionally, an interesting observation is that the EM algorithm takes the longest time when the sample size is small $(N = 200)$ . This suggests that the EM algorithm converges more slowly with smaller sample sizes.

Table 10 The computational time (in seconds) for the VBEM-M, VB, MCMC-R2jags and EM-GDINA algorithms based on LCDM with the $\sigma =0.3$ condition in simulation study 3

5.1 Empirical Example 1

In this example, a fraction subtraction test dataset (de la Torre & Douglas, Reference de la Torre and Douglas2004; Tatsuoka, Reference Tatsuoka, Frederiksen, Glaser, Lesgold and Shafto1990, Reference Tatsuoka2002) was investigated using the DINA model. The VBEM-M algorithm, VB algorithm (implemented in the “variationalDCM” package), MCMC sampling technique (implemented in the “dina” package), and EM algorithm (implemented in the “GDINA” package) were used for the parameter estimation of the DINA model. This test involves 2144 middle school students responding to 15 fraction subtraction items, including five measured attributes: subtract basic fractions, reduce and simplify, separate whole from fraction, borrow from whole, and convert whole to fraction; 536 of 2144 students were chosen for this study (Zhang et al., Reference Zhang, Zhang, Lu and Tao2020). The corresponding $\text {Q}$ -matrix, parameter estimates, and SDs are shown in Table 11.

Table 11 The $\text {Q}$ -matrix and the estimation results of the parameters $\eta $ and $\lambda $ using the VBEM-M algorithm in the empirical example 1

Note: The values outside the parentheses represent the posterior means of the parameters, while the values inside the parentheses indicate the standard deviation.

To facilitate the following item analysis, we transformed the estimates of the intercept and interaction parameters into the traditional estimates of slipping and guessing parameters, as shown in Table 11. Additionally, the comparison of the parameter estimates among the four algorithms can be found in the Supplementary Material. Based on Table 11, we found that the estimates of the five items with the lowest slipping are items 3, 8, 9, 10, and 7, in that order. The estimated values of the slipping parameters for the five items are 0.0395, 0.0480, 0.0652, 0.0664, and 0.0773, respectively. This demonstrates that the five items are less likely to slip than the other ten items. Furthermore, the five items with the highest guessing are items 2, 10, 8, 5, and 13, in that order. For these five items, the estimated guessing parameters are 0.2035, 0.1658, 0.1417, 0.1307, and 0.1293, respectively. Moreover, items 3, 8, and 10 have low slipping parameters and high guessing parameters, indicating that these items are more likely to be correctly guessed. It is worth noting that there is an interesting observation regarding the results for item 1: since $g_1$ is very small and $s_1$ is very large, it is difficult for students who do not master the first attribute to get a correct response by guessing (the probability of a correct response is lower than 0.0200), and even if they do master the first attribute, the probability of a correct response is still only about 0.7000 due to the possibility of slipping.

Based on the results in Table S11 in the Supplementary Material, we investigated the relationship between the VBEM-M algorithm and the other three algorithms in parameter estimation by analyzing the correlations of parameters s and g across these algorithms. The correlations between s estimates from the VBEM-M and VB algorithms are 0.9984, between VBEM-M and MCMC algorithms is 0.9979, and between VBEM-M and EM algorithms is 0.9989. The correlations between the estimators of g calculated using the VBEM-M algorithm and those obtained from the VB, MCMC, and EM algorithms are 0.9488, 0.9552, and 0.8632, respectively. These findings suggest that the VBEM-M algorithm’s parameter estimates align more closely with those from VB and MCMC algorithms, as indicated by the high correlations. In addition, the estimators of the mixing proportions of attribute-mastery patterns, π^ _l for l = 1, ⋯ , 2⁵ = 32, are presented in Figure S2 in the Supplementary Material. Notably, these estimates are highly consistent across the VBEM-M algorithm, VB algorithm, MCMC sampling technique, and EM algorithm. A total of 67% of the examinees were classified into the following four attribute profiles: (1,1,1,0,0), (1,1,1,1,0), (1,1,1,0,1), and (1,1,1,1,1). This suggests that a majority of students have mastered the first three attributes. The computation time for the VBEM-M, VB, MCMC, and EM algorithms were 0.1651 s, 0.1661 s, 11.3820 s, and 0.2870 s, respectively.

5.2 Empirical Example 2

In this section, we analyze the Examination for the Certificate of Proficiency in English (ECPE) dataset based on the LCDM. The ECPE has been widely used in previous research based on the LCDM (e.g., Liu & Johnson, Reference Liu, Johnson, von Davier and Lee2019; Templin & Bradshaw, Reference Templin and Bradshaw2014; Templin & Hoffman, Reference Templin and Hoffman2013; von Davier, Reference von Davier2014b), and it includes 0-1 response data from 2,922 examinees on 28 items. Three attributes are measured: morphosyntactic rules, cohesive rules, and lexical rules. Nine of the 28 items measure two attributes, and the others measure one. The VBEM-M algorithm, VB algorithm (implemented in the “variationalDCM” package), MCMC algotithm (implemented in the “R2jags” package), and EM algorithm (implemented in the “GDINA” package) were used for the parameter estimation of the LCDM model. However, due to space limitations, we only present the estimation results of the VBEM-M method in Tables 12 and 13. The results of the other algorithms can be found in the Supplementary Material.

Table 12 The estimation results of the parameters $\boldsymbol {\eta }$ and $\boldsymbol {\lambda } $ using the VBEM-M algorithm in the empirical example 2

Note: The values outside the parentheses represent the posterior means of the parameters, while the values inside the parentheses indicate the standard deviation.

Table 13 The estimation results of the class membership parameters $\boldsymbol {\pi }$ using the VBEM-M algorithm in the empirical example 2

The outcomes of the VBEM-M algorithm were more similar to those of the VB algorithm and the MCMC algorithm. Please refer to the Supplementary Material for more details. From Table 12, we found that the estimates of the interaction terms are relatively smaller compared to the main effects, indicating that the main effects have a greater influence on the probability of a correct response. Additionally, most of the interaction effects are positive, suggesting that the interactions between skills are more likely to positively affect the probability of a correct response. Furthermore, from the estimates of $\boldsymbol {\pi }$ in Table 13, we can observe that the most prevalent attribute mastery patterns are (0, 0, 0), (0, 0, 1), (0, 1, 1), and (1, 1, 1). This suggests a possible linear hierarchy structure among the skills. Specifically, mastering lexical rules requires mastering cohesive rules first, and mastering morphosyntactic rules is a prerequisite for mastering cohesive rules. This finding is consistent with previous research conclusions (Gierl, Cui, et al., Reference Gierl, Cui and Hunka2007; Gierl, Leighton, et al., Reference Gierl, Leighton, Hunka, Leighton and Gierl2007).