2.1 DINA model

To introduce DINA model formulation, let $i~(i=1,\ldots ,N)$ denote examinees, $j~(j=1,\ldots ,J)$ denote items, and $k~(k=1,\ldots ,K)$ denote attributes. Whether examinee i masters the attribute k is denoted as a binary variable $\alpha _{ik}$ , where $\alpha _{ik}=1$ if examinee i has mastered the attribute k and $\alpha _{ik}=0$ otherwise. We use a K-dimensional binary vector, $\boldsymbol {\alpha }_{i}=(\alpha _{i1},\ldots ,\alpha _{ik},\ldots ,\alpha _{iK})^{\text {T}}$ , to represent the latent attribute profile of examinee i. Therefore, there are $C=2^K$ possible attribute profiles. Let class c denote the bijection index of the attribute profile, where $c=0,\ldots,C-1$ .

Furthermore, a $J \times K $ -dimensional $\text {Q}$ -matrix, $\text {Q}=(\boldsymbol {q}_1,\ldots ,\boldsymbol {q}_j,\ldots ,\boldsymbol {q}_J)^{\text {T}}$ , is introduced to describe which attributes are measured by each item. Here, $\boldsymbol {q}_j^{\text {T}}=(q_{j1},\ldots ,q_{jk},\ldots ,q_{jK})$ is the j-th row of $\text {Q}$ -matrix, corresponding to measured attributes by item j with $q_{ik}=1$ if the attribute k is required by item j, and $q_{ik}=0$ otherwise. Let $Y_{ij}$ denote the binary item response variable for the examinee i to the item j, where $Y_{ij}$ equals 1 if the ith examinee correctly answered the item j, and 0 otherwise.

The DINA model assumes an examinee need to master all the required attributes for item j to have a high chance to answer that item correctly. To reflect this assumption, the ideal response variable $\eta _{ij}=\prod _{k=1}^{K}\alpha _{ik}^{q_{jk}}$ is defined to formulate the item response model. Here, $\eta _{ij}=1$ indicates that examinee i possesses all the required attributes for item j; otherwise, $\eta _{ij}=0$ . In addition, two item parameters are introduced for each item j: the slipping parameter $s_j$ and the guessing parameter $g_j$ which can be formally defined by

(1) $$ \begin{align} s_j&=P(Y_{ij}=0|\eta_{ij}=1),\ g_j=P(Y_{ij}=1|\eta_{ij}=0). \end{align} $$

Then, the probability of a correct response for the DINA can be expressed as

(2) $$ \begin{align} P(Y_{ij}=1|\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_{c},g_j,s_j)=g_j^{1-\eta_{cj}}(1-s_j)^{\eta_{cj}}, \end{align} $$

where $\boldsymbol {\alpha }_{c}$ denotes one possible attribute profile value and $\eta _{cj}$ is the corresponding ideal response variable.

2.2 Attribute hierarchy and $\text {Q}$ -matrix

One type of latent structure with CDM is the attribute hierarchy, which may be commonly considered when defining attributes within many educational and psychological tests. Based on the nature of assessment, attributes may have different types of dependency, where one attribute may be a prerequisite for another attribute. Four types of popular hierarchical structure are discussed by Leighton et al. (Reference Leighton, Gierl and Hunka2004), which are presented as four directed acyclic graphs (DAGs) in Figure 1. For a DAG, each node represents an attribute, and a directed edge between two nodes indicates a prerequisite relationship, where the attribute represented by the source node is required for mastering the attribute represented by the target node. For instance, with the linear hierarchy represented by Figure 1(a), the directed edge $1\rightarrow 2$ indicates that the mastery of $\alpha _1$ is a prerequisite for mastering $\alpha _2$ .

Figure 1 Four attribute hierarchy structures, arranged from left to right: (a) $\text {G}_1$ Linear; (b) $\text {G}_2$ Convergent; (c) $\text {G}_3$ Divergent; and (d) $\text {G}_4$ Unstructured.

In this article, we adopt the same notation for describing the attribute hierarchy as Chen and Wang (Reference Chen and Wang2023). Let $ V $ and E denote the set of attributes and the directed edges, respectively, and $G(V;E)$ denote the attribute hierarchy structure. Mathematically, one $G(V;E)$ can be presented by its corresponding adjacency matrix, $\text {G}$ . For example, let $G(V;E)$ represent the linear hierarchy in Figure 1(a). With this structure, $V=\{1,2,3,4\}$ and $E=\{1 \rightarrow 2,2 \rightarrow 3, 3 \rightarrow 4\}$ . This particular $G(V;E)$ can be denoted by its adjacency matrix $\text {G}$ as

$$ \begin{align*}\text{G}=\begin{bmatrix} 0 & 1 & 0 &0\\ 0& 0 & 1&0\\ 0& 0& 0&1\\ 0& 0& 0&0\\ \end{bmatrix}, \end{align*} $$

where $\text {G}(k,k^{'})=1$ represents there is a directed edge $k \rightarrow k'$ and 0 otherwise. When there is no hierarchical structure between any attributes, all elements in $\text {G}$ are 0, and we denote this structure as $\text {G}_{null}$ . Because of this one-to-one correspondence between a DAG and its adjacency matrix $\text {G}$ , we will just use the adjacency matrix $\text {G}$ to denote its corresponding hierarchical structure $G(V,E)$ .

In addition to the direct relationships discussed above, there are also indirect dependencies among attributes. To better capture both direct and indirect dependencies between attributes, the reachability matrix $\text {R}$ is introduced, where $\text {R}(k,k^{'})=1$ represents that there exists a path from k to $k'$ , and it can be obtained through Boolean addition and multiplication of the adjacency matrix (Tatsuoka, Reference Tatsuoka2012). The reachability matrix $\text {R}$ of Figure 1(a) can be expressed as

$$ \begin{align*}\text{R}=\begin{bmatrix} 1 & 1 & 1 &1\\ 0& 1 & 1&1\\ 0& 0& 1&1\\ 0& 0& 0&1\\ \end{bmatrix}. \end{align*} $$

Note that different hierarchical structures may lead to the same reachability matrix. For example, consider the structure $ \text {G}^* $ , corresponding to $ V^* = \{1, 2, 3, 4\} $ and $E^* = \{1 \rightarrow 2, 2 \rightarrow 3, 3 \rightarrow 4, 1 \rightarrow 4\}$ . The reachability matrix of $ \text {G}^* $ will be the same as that of $ \text {G} $ for Figure 1(a), meaning that both exhibit the same dependencies between attributes. However, it is clear that the structure of $ \text {G} $ in Figure 1(a) is simpler. Therefore, our proposed Bayesian estimation method estimates the transitive reduction of $\text {G}$ , which refers to a subgraph of $ \text {G} $ with the fewest edges that maintains the same reachability as $\text {G}$ (Chen & Wang, Reference Chen and Wang2023).

2.2.2 $\text {Q}$ -matrix under $\text {G}$

Another important latent structure under CDM is the $\text {Q}$ -matrix. When attribute hierarchy exists, we assume that the hierarchical structure also affects the space of the $\text {Q}$ -matrix. In other words, the $\boldsymbol {q}$ -vector needs to reflect the underlying attribute hierarchy as well. For example, considering the linear structure $\text {G}_1$ , we need to make the notations consistent in Figure 1(a), $\mathbf {q} = (0, 1,0,0)$ equivalent to $\mathbf {q} = (1, 1,0,0)$ , indicating that an item measuring $\alpha _2$ also needs to measure $\alpha _1$ under this linear structure. Therefore, we introduce reduced $\text {Q}$ -matrix $\text {Q}_r$ (Tatsuoka, Reference Tatsuoka2012), which including all possible $\boldsymbol {q}$ -vectors under the hierarchical structure. In this case, the total $2^K-1$ possible $\boldsymbol {q}$ -vectors under $\text {G}_{null}$ are reduced to four possible $\boldsymbol {q}$ -vectors under the linear structure $\text {G}_1$ . In this study, a structured $\text {Q}$ -matrix constructed from the possible $\boldsymbol {q}$ -vectors in $\text {Q}_r(\text {G})$ is used to reflect the attribute hierarchical structure. The above relationship is shown by Equation (3), where each column in a $\text {Q}$ -matrix represents a possible $\boldsymbol {q}$ -vector:

(3) $$ \begin{align} \text{Q}_{all}^{\text{T}} = \left[\begin{array}{lllllllllllllll} 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & 1 & 1 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \end{array}\right] \stackrel{reduced}{\longrightarrow} \left(\begin{array}{lllll} 1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{array}\right) \stackrel{\triangle}{=}\text{Q}_r (\text{G}_1)^{\text{T}}. \end{align} $$

3.1 Bayesian formulation

To start with, denote item response matrix for N examinees across J items as $\textbf {Y}=(\boldsymbol {Y}_{1},\ldots ,\boldsymbol {Y}_{i},\ldots ,\boldsymbol {Y}_{N})^{\text {T}}$ , where $\boldsymbol {Y}_{i}=(y_{i1},\ldots ,y_{ij},\ldots ,y_{iJ})^{\text {T}}$ , $i=1,\ldots ,N$ , is the response vector of examinee i to J items. In addition, let $\textbf {A}=(\boldsymbol {\alpha }_{1},\ldots ,\boldsymbol {\alpha }_{i},\ldots ,\boldsymbol {\alpha }_{N})^{\text {T}}$ be the $N\times K$ matrix of attribute profiles of all examinees. Denote the parameter $\boldsymbol {\pi }=(\pi _1,\ldots ,\pi _c,\ldots ,\pi _{C})$ as the latent class probability vector, that is, $\pi _c=P(\boldsymbol {\alpha }_i=\boldsymbol {\alpha }_c)$ , satisfying $\sum \pi _c=1$ . Let $\boldsymbol {s}=(s_1,\dots ,s_J)^{T}$ and $\boldsymbol {g}=(g_1,\dots ,g_J)^{T}$ denote the slipping and guessing parameters for all items.

Given the conditional independent assumption, the joint distribution of $\mathbf {Y}$ given parameters $\boldsymbol {s},\boldsymbol {g},\boldsymbol {\pi }, \boldsymbol {A}, \text {Q}$ , and $\text {G}$ can be written as

(4) $$ \begin{align} p(\mathbf{Y}|\boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G})=\prod_{i=1}^N\prod_{j=1}^J\sum_{c=1}^C \boldsymbol{\pi}_c ~p(Y_{ij}|\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_c,g_j,s_j, \boldsymbol{q}_j,\text{G}), \end{align} $$

where

(5) $$ \begin{align} \begin{aligned} p(Y_{ij}|\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_c,g_j,s_j,\boldsymbol{q}_j, \text{G})&=\left(g_j^{1-\eta_{cj}}(1-s_j)^{\eta_{cj}}\right)^{y_{ij}} \left((1-g_j)^{1-\eta_{cj}}s_j^{\eta_{cj}}\right)^{1-y_{ij}},\\ \eta_{cj}&=\prod_{k=1}^{K}\alpha_{ck}^{q_{jk}}. \end{aligned} \end{align} $$

To estimate the DINA model item parameters $\boldsymbol {s}$ and $\boldsymbol {g}$ , the examinee parameters $\boldsymbol {A}$ and $\boldsymbol {\pi }$ , and the two latent structures $\text {G}$ and $\text {Q}$ , we first outline the joint posterior distribution of all model parameters.

Let $p(s_j,g_j)$ , $p(\boldsymbol {\alpha }_i\mid \boldsymbol {\pi },\text {G})$ , $p(\boldsymbol {\pi }\mid \text {G})$ , $~p(\text {Q}\mid \text {G})$ , and $p(\text {G})$ be the corresponding prior of $(s_j,g_j)$ , $\boldsymbol {\alpha }_i$ , $\boldsymbol {\pi }$ , $\text {Q}$ , and $\text {G}$ , respectively. The joint posterior distribution of all model parameters is

(6) $$ \begin{align} \begin{aligned} &p(\boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G}\mid \mathbf{Y}) \\&\quad\propto ~p(\mathbf{Y}\mid \boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G}) p(\boldsymbol{s},\boldsymbol{g})~p(\mathbf{A}\mid \boldsymbol{\pi},\text{G})~p(\boldsymbol{\pi}\mid \text{G})~p(\text{Q}\mid \text{G})~p(\text{G})\\ &\quad\propto ~p(\mathbf{Y}\mid \boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G})\left(\prod_{j=1}^J p(s_j,g_j)\right) \left(\prod_{i=1}^N p(\boldsymbol{\alpha}_i \mid \boldsymbol{\pi},\text{G})\right)~ p(\boldsymbol{\pi}\mid \text{G})~p(\text{Q}\mid \text{G})~p(\text{G}), \end{aligned} \end{align} $$

where $p(\mathbf {Y}\mid \boldsymbol {s},\boldsymbol {g},\boldsymbol {\pi },\boldsymbol {A},\text {Q},\text {G})$ is the joint distribution function (4). In the following section, we specify the prior distributions of each parameter, followed by the inference of their posterior distributions.

3.2 Posterior inference

Firstly, the prior of $s_j$ and $g_j$ follows the Beta distribution with parameter $a_{s0}$ , $b_{s0}$ , $a_{g0}$ , and $b_{g0}$ :

(7) $$ \begin{align} p(s_j,g_j)\propto s_j^{a_{s0}-1} (1-s_j)^{b_{s0}-1}g_j^{a_{g0}-1} (1-g_j)^{b_{g0}-1} \mathcal{I}(0<s_j+g_j<1). \end{align} $$

Then we can derive that posterior distribution of $s_j$ is a Beta distribution with parameter $a_{sj}$ and $b_{sj}$ as following:

(8) $$ \begin{align} s_j \mid \textbf{Y},\boldsymbol{g},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G} \sim \operatorname{Beta}(a_{sj},b_{sj}) \mathcal{I}(0<s_j<1-g_j), \end{align} $$

where

$$ \begin{align*} \begin{aligned} a_{sj} = \sum_{i=1}^{N} \eta_{ij}(1-y_{ij}),\quad b_{sj} = \sum_{i=1}^{N} \eta_{ij}y_{ij}. \end{aligned} \end{align*} $$

And the posterior distribution of $g_j$ is a Beta distribution with parameter $a_{gj}$ and $b_{gj}$ as following:

(9) $$ \begin{align} g_j \mid \textbf{Y},\boldsymbol{s},\boldsymbol{\pi},\boldsymbol{A},\text{Q},\text{G} \sim \operatorname{Beta}(a_{gj},b_{gj}) \mathcal{I}(0<g_j<1-s_j), \end{align} $$

where

$$ \begin{align*} \begin{aligned} a_{gj} = \sum_{i=1}^{N} (1-\eta_{ij})y_{ij},\quad b_{gj} = \sum_{i=1}^{N} (1-\eta_{ij})(1-y_{ij}). \end{aligned} \end{align*} $$

Following Culpepper (Reference Culpepper2015), we assign Beta(1,1) priors to both $s_j$ and $g_j$ in this article, that is, $a_{s0} = b_{s0} = a_{g0} = b_{g0} = 1$ .

Secondly, the prior of $\boldsymbol {\pi }$ under structure $\text {G}_{null}$ is a Dirichlet distribution with parameter $\boldsymbol {\delta }_0$ :

(10) $$ \begin{align} p(\boldsymbol{\pi} \mid \boldsymbol{\delta}_0) \propto \prod_{c=0}^{C-1} \pi_{c}^{ \delta_{0c}}, 0 \leq \pi_c\leq 1, \sum \pi_c=1. \end{align} $$

However, as discussed in Section 2.2, the existence of a hierarchical structure among attributes imposes constraints on the attribute profiles. Specifically, some attribute profiles become impermissible under the hierarchy structure $\text {G}$ . This restriction reduces the set of permissible attribute profiles, thereby influencing the distribution of $\boldsymbol {\pi }$ . Moreover, since $\text {G}$ is allowed to change during the estimation process, the distributions of $\boldsymbol {\alpha }_i$ also change accordingly. Therefore, when incorporating a hierarchical structure $\text {G}$ , it is necessary to adjust the corresponding priors to reflect the structural constraints imposed by $\text {G}$ .

Let $L_G$ denote the set of permissible attribute profiles under a specific hierarchical structure $\text {G}$ and $C_G=|L_G|$ denote the number of permissible attribute profiles. Define $\boldsymbol {\pi }_G$ as the current latent class probability vector. Thus, the prior of $\boldsymbol {\pi }_G$ now follows an aggregated Dirichlet distribution with parameters $\boldsymbol {\delta }_0^G$ :

(11) $$ \begin{align} p(\boldsymbol{\pi}_G \mid \boldsymbol{\delta}_0^G) \propto \prod_{\boldsymbol{\alpha}_{\tilde{c}} \in L_G} \pi_{\tilde{c}}^{ \delta^G_{0\tilde{c}}}, \end{align} $$

where

$$ \begin{align*} \pi_{\tilde{c}}=\sum_{\boldsymbol{\alpha}_c \stackrel{G}{=} \boldsymbol{\alpha}_{\tilde{c}}; c \in C} \pi_c, \quad \delta_{0\tilde{c}}^G=\sum_{\boldsymbol{\alpha}_c \stackrel{G}{=} \boldsymbol{\alpha}_{\tilde{c}}; c \in C} \delta_c. \end{align*} $$

Based on the prior in Eq. (11) and the joint posterior distribution in Eq. (6), we can derive that posterior distribution of $\boldsymbol \pi _G$ is a Dirichlet distribution with parameter $\boldsymbol {\delta }_G$ as follows:

(12) $$ \begin{align} \boldsymbol{\pi}_G \mid \textbf{Y},\boldsymbol{s},\boldsymbol{g}, \textbf{A},\text{Q},\text{G} \sim \operatorname{Dirichlet}(\boldsymbol{\delta}_G), \end{align} $$

where

$$ \begin{align*} \delta^G_{\tilde{c}}= \sum_{i=1}^{N}\mathcal{I}(\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_{\tilde{c}})+\delta_{0\tilde{c}}^G. \end{align*} $$

We set $\boldsymbol {\delta }_0 = (1, \ldots , 1)_L$ in our study, following the setting in Culpepper (Reference Culpepper2015).

Thirdly, the prior of $\boldsymbol {\alpha }_i$ follows a categorical distribution with parameter $\boldsymbol {\pi }$ under structure $\text {G}_{null}$ :

(13) $$ \begin{align} p\left(\boldsymbol{\alpha}_i | \boldsymbol{\pi}\right) = \prod \pi_{c}^{\mathcal{I}(\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_{c})}, c=0,\ldots, C-1. \end{align} $$

Similar to parameter $\boldsymbol {\pi }$ , the prior of $\boldsymbol {\alpha }_i$ also needs to be adjusted according to the hierarchy structure $\text {G}$ . Accordingly, the prior of $\boldsymbol {\alpha }_i$ follows a collapsed categorical distribution with parameter $\boldsymbol {\pi }_G$ under structure G:

(14) $$ \begin{align} p\left(\boldsymbol{\alpha}_i | \boldsymbol{\pi}_G\right) = \prod_{\tilde{c}\in L_G} \pi_{\tilde{c}}^{\mathcal{I}(\boldsymbol{\alpha}_i=\boldsymbol{\alpha}_{\tilde{c}})}. \end{align} $$

Based on the prior in Eq. (14) and the joint posterior distribution in Eq. (6), we can derive that posterior distribution of $\boldsymbol {\alpha }_i$ is a categorical distribution with parameter $\boldsymbol {\rho }_{i}^G=(\rho _{i0}^G,\ldots ,\rho _{iC_G-1}^G)$ as follows:

(15) $$ \begin{align} \boldsymbol{\alpha}_i\mid \textbf{Y},\boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi}_G,\text{Q},\text{G} \sim \operatorname{Categorical}(\boldsymbol{\rho}_{i}^G), \end{align} $$

where

$$ \begin{align*} \rho_{i\tilde{c}}^G \propto \pi_{\tilde{c}}^G \prod_{j=1}^J \left(g_j^{1-\eta_{\tilde{c}j}}(1-s_j)^{\eta_{\tilde{c}j}}\right)^{y_{ij}} \left((1-g_j)^{1-\eta_{\tilde{c}j}}s_j^{\eta_{\tilde{c}j}}\right)^{1-y_{ij}}. \end{align*} $$

Fourthly, as for the prior of the $\text {Q}$ -matrix, we adopted the approach in Chung (Reference Chung2019) by introducing additional parameters for each Q-matrix entry and made corresponding adjustments. With K attributes, there are $H=2^K-1$ possible $\boldsymbol {q}$ -vector patterns for each item under structure $\text {G}_{null}$ , which is denoted as $\boldsymbol {\epsilon }_{H \times K}=\left (\epsilon _{h k}\right )_{H \times K}$ . Suppose $\boldsymbol {q}_j$ following a categorical distribution with parameter $\boldsymbol {\theta }_{0j}=(\theta _{0j1},\ldots ,\theta _{0jH})$ :

(16) $$ \begin{align} p(\boldsymbol{q}_j \mid \boldsymbol{\theta}_{0j}) \propto \prod_{h=1}^H \theta_{0jh}^{\mathcal{I}(\boldsymbol{q}_j=\boldsymbol{\epsilon}_{h})}, \end{align} $$

$$ \begin{align*} \boldsymbol{\theta}_{0j} = g(\boldsymbol{\Phi}_{j}),~~ \boldsymbol{\Phi}_{j} =(\boldsymbol{\phi}_{j1},\ldots,\boldsymbol{\phi}_{jH})^{\text{T}}, \end{align*} $$

where $\boldsymbol {\phi }_{jh} =(\phi _{jh1},\ldots ,\phi _{jhK}), h\in (1,\ldots ,H)$ . We define $\phi _{jhk}=P(q_{jk}=1)$ and the prior for $\phi _{jhk}$ is

$$ \begin{align*} \phi_{jhk} \sim \operatorname{Beta}(1,1). \end{align*} $$

Therefore, the conditional posterior for each element $\phi _{jhk}$ is distributed as

(17) $$ \begin{align} \phi_{jh k}\mid q_{jk} \sim \operatorname{Beta}\left(1+q_{j k}, 2-q_{j k} \right). \end{align} $$

Moreover, the prior for $\boldsymbol {\theta }_{0j}$ is then modeled as

(18) $$ \begin{align} \begin{aligned} \boldsymbol{\theta}_{0j}=g( \boldsymbol{\Phi}_j)& =(g(\boldsymbol{\phi}_{j1}),\ldots,g(\boldsymbol{\phi}_{jH})\propto \left(\prod_{k=1}^K \phi_{j1 k}^{\epsilon_{1 k}}\left(1-\phi_{j1 k}\right)^{1-\epsilon_{1 k}},\ldots, \prod_{k=1}^K \phi_{jh k}^{\epsilon_{h k}}\right. \\ &\left.\quad \times\left(1-\phi_{jh k}\right)^{1-\epsilon_{h k}}, \ldots, \prod_{k=1}^K \phi_{jH k}^{\epsilon_{H k}}\left(1-\phi_{jH k}\right)^{1-\epsilon_{H k}}\right). \end{aligned} \end{align} $$

Considering that in the DINA model, restricting $\boldsymbol {\alpha }$ is equivalent to restricting the $ \text {Q} $ -matrix, we do not impose any additional structural constraints on the $ \text {Q} $ -matrix during its update process. Therefore, based on the prior in Eq. (18) and the joint posterior distribution in Eq. (6), we can derive that posterior distribution of $\boldsymbol {q}_j$ is a categorical distribution with parameter $\boldsymbol {\theta }_j=(\theta _{j1},\ldots ,\theta _{jH})$ as follows:

(19) $$ \begin{align} \boldsymbol{q}_j\mid \textbf{Y},\boldsymbol{s},\boldsymbol{g},\boldsymbol{\pi},\text{G} \sim \operatorname{Categorical}(\boldsymbol{\theta}_{j}), \end{align} $$

where

(20) $$ \begin{align} \theta_{jh}=\frac{\gamma_{jh}}{\sum_{h=1}^H\gamma_{jh}}, \end{align} $$

(21) $$ \begin{align} \gamma_{jh} = p(\boldsymbol{\phi}_{jh}) \prod_{i=1}^N \left(g_j^{1-\eta_{ij}}(1-s_j)^{\eta_{ij}}\right)^{y_{ij}} \left((1-g_j)^{1-\eta_{ij}}s_j^{\eta_{ij}}\right)^{1-y_{ij}}. \end{align} $$

Finally, for the hierarchical structure $\text {G}$ , we adopt a uniform prior, which means for any two possible transitive reduction structure $\text {G}$ and $\text {G}^*$ , we have $p(\text {G})=p(\text {G}^*)$ . As we are not able to get the closed form of the posterior distribution of G, similar to Chen and Wang (Reference Chen and Wang2023), the MH sampling method is utilized to update $\text {G}$ . First, we need to sample a new structure $\text {G}^*$ . Note that the new $\text {G}^*$ can only be obtained by either removing an existing edge or adding an edge to the current structure $\text {G}$ . Denote p as the probability of adding an edge, indicating $1-p$ as the probability of removing an edge. We first use the probability p to determine whether the new structure $\text {G}^*$ involves adding or removing an edge, and then randomly generate a possible $\text {G}^*$ . It is worth to note that this sampling procedure only exists between transitive reduction $\text {G}$ s, meaning that each time a new $\text {G}^*$ is sampled, it must have a different reachability matrix from $\text {G}$ . According to this rule, each time an edge is added, it must be selected from the edges between two attributes that do not have any direct or indirect relationship. Such a pair of attributes can be found through $0$ s in the off-diagonal entries in the reachability matrix. Let $T\left (\text {G}, \text {G}^*\right )$ denote the transition probability from $\text {G}$ to the proposed new $\text {G}^*$ . Define

(22) $$ \begin{align} T\left(\text{G}, \text{G}^*\right)=\frac{1}{2 \times \text { number of pairs of non-connected nodes }}, \end{align} $$

when adding an edge, and

(23) $$ \begin{align} T\left(\text{G}, \text{G}^*\right)=\frac{1}{\text { number of existing edges }}, \end{align} $$

when removing an edge. Then the accept ratio of $\text {G}^*$ is

(24) $$ \begin{align} \begin{aligned} r & =\min \left\{1, \frac{p\left(\text{G}^*\mid \textbf{Y}, \boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \mathbf{Q}\right) T\left(\text{G}^*, \text{G}\right)}{p\left(\text{G}\mid \textbf{Y}, \boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \mathbf{Q}\right) T\left(\text{G}, \text{G}^*\right)}\right\} \\ & =\min \left\{1, \frac{p\left(\textbf{Y} \mid \text{G}^*,\boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \text{Q}\right) p\left(\text{G}^*\right) T\left(\text{G}^*, \text{G}\right)}{p(\textbf{Y} \mid \text{G},\boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \text{Q}) p(G) T\left(\text{G}, \text{G}^*\right)}\right\}\\ &=\min \left\{1, \frac{p\left(\textbf{Y} \mid \text{G}^*,\boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \mathbf{Q}\right) T\left(\text{G}^*, \text{G}\right)}{p(\textbf{Y} \mid \text{G}, \boldsymbol{s},\boldsymbol{g}, \mathbf{A}, \text{Q}) T\left(\text{G}, \text{G}^*\right)}\right\}. \end{aligned} \end{align} $$

Then, a random number is drawn from a uniform distribution and compared to the acceptance ratio r to determine whether to accept the proposed $\text {G}^*$ .

Our proposed MH within Gibbs sampling procedure is summarized in Table 1, and T represent the length of the Markov chain. The prespecified probability p can be adjusted to improve the acceptance ratio of the newly proposed structure $\text {G}^*$ . In this article, we set $p=0.5$ .

Table 1 MH-Gibbs sampling procedure

4.1 Mini-batch method in sampling procedure

We propose to use the mini-batch method, where, during each iteration of the sampling algorithm (Table 1), we randomly select a subset of samples to update $\boldsymbol {s}$ , $\boldsymbol {g}$ , $\boldsymbol {Q}$ , and $\text {G}$ . Note that $\boldsymbol {\alpha }$ and $\boldsymbol {\pi }$ cannot be updated using sub-samples, as they are parameters specific to each individual sample. The mini-batch method is a widely used technique in machine learning, commonly applied in stochastic gradient descent. The method works by updating model parameters using only a subset of the samples during each iteration. In recent years, the mini-batch method has also found applications in the Bayesian field, such as in MH (Wu et al., Reference Wu, Rachel Wang and Wong2022), Gibbs sampling (De Sa et al., Reference De Sa, Chen and Wong2018), Hamiltonian Monte Carlo (HMC) (Zou & Gu, Reference Zou and Gu2021), and variational Bayes (Hoffman et al., Reference Hoffman, Blei, Wang and Paisley2013). This approach offers two main advantages. First, it improves the computational efficiency of the algorithm. Second, the mini-batch method introduces some noise, thereby increasing the diversity of the sampling and helping to prevent rapid convergence to local optima.

When implementing mini-batch with our proposed algorithm, one needs to determine the mini-batch size. Common practice in machine learning involves selecting batch sizes that are powers of two, such as 32, 64, 128, or 256, due to their compatibility with CPU and GPU memory architectures (Keskar et al., Reference Keskar, Mudigere, Nocedal, Smelyanskiy and Tang2016). Keskar et al. (Reference Keskar, Mudigere, Nocedal, Smelyanskiy and Tang2016) have discussed that large-batch sample size tend to converge to sharp minimizers, while small-batch sample size consistently converge to flatter minimizers. This suggests that using a larger batch size does not necessarily lead to better estimation performance. In contrast, a relatively smaller batch size may introduce more stochasticity, which can help the algorithm escape from local optima. Given these considerations, we will specifically explore the impact of different mini-batch size to the model estimation results in the simulation section.

4.2 Label switch

Before computing the final estimations from the sampling procedure, it is necessary to address the issue of label switching. In the context of Bayesian methods, label switching is a common concern that may arise both within a single Markov chain and across multiple chains. In this article, the unknown Q-matrix may lead to label switching issues among several attributes, which also impacts $\boldsymbol {\alpha }$ and $\text {G}$ . Therefore, it is crucial to impose a relabeling process on the iterative data. Denote the $J \times K \times T$ -dimension $\mathcal {Q}$ as the iteration samples of Q. Following the relabeling method in Erosheva and Curtis (Reference Erosheva and Curtis2017), we define $\overline {\mathcal {Q}}$ as the average of the sampled $\mathcal {Q}$ , which serves as the first reference. For each $\text {Q}^{(t)} \in \mathcal {Q}$ , we leverage the Hungarian algorithm to address the relabeling issue. Specifically, we construct a $K \times K$ cost matrix $C^{(t)} = \{c_{kk'}^{(t)}\}$ , where

$$\begin{align*}c_{kk'}^{(t)} = \big\| \mathbf{q}_{k'}^{(t)} - \overline{\mathbf{q}_k} \big\|_1, \quad k,k'=1,\ldots,K,\; t=1,\ldots,T. \end{align*}$$

Let $\mathcal {P}$ denote the set of all permutations of $\{1,\ldots ,K\}$ . Then, the relabeled $\text {Q}^{(t)}$ , denoted as $\text {Q}_R^{(t)}$ , is obtained by solving the following assignment problem:

(25) $$ \begin{align} p^{\ast} = \arg\min_{p \in \mathcal{P}} \sum_{k=1}^K c_{k,p(k)}^{(t)}, \end{align} $$

(26) $$ \begin{align} \text{Q}_R^{(t)} = \text{Q}_{p^{\ast}}^{(t)}, \quad t=1,\ldots,T. \end{align} $$

Note that the above steps need to be repeated. In other words, after obtaining the new relabeled $\mathcal {Q}_R$ , we need to calculate the new average matrix $\overline {\mathcal {Q}}_R$ and repeat the relabel procedure until $\mathcal {Q}_R$ converge. In addition, the same permutation must be applied to both $\boldsymbol {\alpha }$ and $\text {G}$ simultaneously, as they are also subject to the same relabel issues involved in Q matrix.

4.3 Final estimation summary

5.1 Design

We consider a set up that a total of $K=4$ attributes were assessed, and we fix the slipping and guessing parameters in a practical level as $s_j=g_j=0.2$ . The summary of the five factors considered and their corresponding level is presented in Table 2. These factors were fully crossed, yielding 40 conditions. We conducted four chains in each replication, with the chain length of 20,000 and the burn in value of 10,000, and each simulation repeated 50 times.

Table 2 Simulation factors and levels

The two test lengths represent relatively short and medium-length assessments. Recognizing that estimating a high-dimensional, unknown Q-matrix necessitates a larger sample size for accuracy, we increased the sample size paired with the increase of test length in our simulation design. This adjustment reflects both practical considerations and the need for a sufficiently large sample size. The four hierarchical structures of the attributes, as well as studied in previous related work (Chen & Wang, Reference Chen and Wang2023), were therefore also considered in our study.

As for the mini-batch size, as noted in Section 3.3, common practice in machine learning involves selecting batch sizes that are powers of two. To balance convergence speed and computational efficiency, we avoided batch sizes that are excessively small or large. Consequently, we chose batch sizes ranging from one-quarter to one-half of the total sample size. Specifically, for $N=500$ with $J=20$ , we selected batch sizes of 128 and 256; for $N=1,000$ with $J=40$ , we chose batch sizes of 128, 256, and 512.

Given the findings in Chen and Wang (Reference Chen and Wang2023), for edges that do not exist and have no indirect relationships between attributes, the estimated values are close to 0. For non-existing edges with indirect relationships between attributes, the estimated values range between 0.2 and 0.4. Therefore, when selecting the cut-off value, we choose 0.2, 0.3, and 0.4.

Finally, following the design in Chen and Wang (Reference Chen and Wang2023), we also considered a balanced and unbalanced attribute profile distribution. For the balanced distribution, the attribute profiles were randomly generated from one of its permissible patterns. For the unbalanced distribution, we first generated $\boldsymbol {\alpha }^*_{i}=(\alpha ^*_{i1},\ldots ,\alpha ^*_{ik},\ldots ,\alpha ^*_{iK})^{\text {T}}$ from a multivariate normal distribution $\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*} $$

the off-diagonal elements of $\Sigma _{K\times K}$ are $\sigma =0.5$ . Then the relationships between the attribute profile $\boldsymbol {\alpha }_i$ and $\boldsymbol {\alpha }^*_i$ can be expressed as $\alpha _{ik}=1$ if $\alpha _{ik}^*>0$ and $\alpha _{ik}=0$ otherwise.

5.2 $\text {Q}$ -matrix design

The establishment of parameter identifiability is a prerequisite for the accurate estimation and reliable interpretation of CDMs. To this end, it is essential to impose a set of constraints on the $\text {Q}$ -matrix. Previous studies have discussed the identifiability conditions of parameters in CDMs (Chen et al., Reference Chen, Liu, Xu and Ying2015; Chen et al., Reference Chen, Culpepper, Chen and Douglas2018; Gu & Xu, Reference Gu and Xu2023; Ma et al., Reference Ma, Ouyang and Xu2023). In particular, Gu and Xu (Reference Gu and Xu2023) explicitly addressed the identifiability conditions when a hierarchical structure exists among attributes and the $\text {Q}$ -matrix is unknown. Therefore, in this study, we followed the identifiability conditions proposed by Gu and Xu (Reference Gu and Xu2023) when generating the $\text {Q}$ -matrix. Specifically, they defined the “densified” version $\mathcal {D}^{\text {G}}(\text {Q})$ and the “sparsified” version $\mathcal {S}^{\text {G}}(\text {Q})$ of the $\text {Q}$ -matrix. The “densified” version $\mathcal {D}^{\text {G}}(\text {Q})$ corresponds to what we referred to earlier as the structured $\text {Q}$ -matrix in Section 2.2, which is defined as: for any $q_{j,k'} = 1$ and $k \to k'$ , set $q_{j,k} = 1$ . Conversely, the “sparsified” version $\mathcal {S}^{\text {G}}(\text {Q})$ is defined as: for any $q_{j,k'} = 1$ and $k \to k'$ , set $q_{j,k} = 0$ . For example, under the structure $\text {G}_1$ , given a q-vector $\boldsymbol {q} = (0,1,1,0)$ , the “densified” version $\mathcal {D}^{\text {G}}(\boldsymbol {q})$ becomes $(1,1,1,0)$ , whereas the “sparsified” version $\mathcal {S}^{\text {G}}(\boldsymbol {q})$ becomes $(0,0,1,0)$ . Moreover, under $G_1$ , these three versions are equivalent, that is,

$$\begin{align*}(0,1,1,0) \;\stackrel{G_1}{\sim}\; (1,1,1,0) \;\stackrel{G_1}{\sim}\; (0,0,1,0). \end{align*}$$

Denote the $\text {Q}$ -matrix contains the following two components:

$$ \begin{align*}\mathbf{Q} = \left( \begin{array}{c} \mathbf{Q}^0_{K\times K} \\ \mathbf{Q}^{\star}_{(J-K)\times K} \end{array} \right). \end{align*} $$

Specifically, we generate the Q-matrix satisfying the following three identifiability conditions: (A) the sub-matrix $\text {Q}^0$ is equivalent to the identity matrix $I_K$ , that is, $\text {Q}^0 \stackrel {\text {G}}{\sim }I_K$ ; (B) the $\mathcal {S}^{\text {G}}(\text {Q})$ has at least three ones in each column; and (C) the $\mathcal {D}^{\text {G}}\left (\text {Q}^{\star }\right )$ contains K distinct column vectors.

Moreover, once we generated $\text {Q}$ -matrices that satisfied the identifiability conditions, we uniformly adopted their corresponding “densified” versions, that is, the structured Q-matrices. It is worth noting that during the 50 replications under each simulation condition, we generated 50 distinct Q-matrices. This design may reduce comparability across different conditions. However, it allows us to better evaluate the performance and robustness of our algorithm across different Q-matrices under each condition.

5.3 Evaluation criteria

We evaluate the proposed algorithm through the following perspectives. First, the potential scale reduction factor (PSRF), denoted as $\hat {R}$ (Gelman & Rubin, Reference Gelman and Rubin1992), was used to evaluate the convergence of the Gibbs algorithm across different simulation conditions. When the $\hat {R}$ is less than 1.1, the Markov chain is considered to have reached convergence. For each replication, four Markov chains were run. The chain that achieved the lowest deviance information criterion (DIC; Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and Van Der Linde2002), a widely used model selection criterion in the Bayesian framework, was considered the best and was selected for subsequent analysis.

Second, we used the average bias and root mean squared error (RMSE) to evaluate the accuracy of the estimation for the item parameters $\boldsymbol {s}$ , $\boldsymbol {g}$ , and attribute profile membership parameter $\boldsymbol {\pi }$ . These evaluation indexes are defined as

$$ \begin{align*} & \operatorname{ABias}(\boldsymbol{s})=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{R} \sum_{r=1}^R\left({\hat{s_{j }}}^{(r)}-s_{j }\right),\quad \operatorname{ARMSE}(\boldsymbol{s})=\frac{1}{J}\sum_{j=1}^{J}\sqrt{\frac{1}{R} \sum_{r=1}^R\left({\hat{s_{j }}}^{(r)}-s_{j }\right)^2},\\ & \operatorname{ABias}(\boldsymbol{g})=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{R} \sum_{r=1}^R\left({\hat{g_{j }}}^{(r)}-g_{j }\right),\quad \operatorname{ARMSE}(\boldsymbol{g})=\frac{1}{J}\sum_{j=1}^{J}\sqrt{\frac{1}{R} \sum_{r=1}^R\left({\hat{g_{j }}}^{(r)}-g_{j }\right)^2},\\ & \operatorname{ABias}(\boldsymbol{\pi})=\frac{1}{C}\sum_{c=1}^{C}\frac{1}{R} \sum_{r=1}^R\left({\hat{\pi_{c }}}^{(r)}-\pi_{c }\right),\quad \operatorname{ARMSE}(\boldsymbol{\pi})=\frac{1}{C}\sum_{c=1}^{C}\sqrt{\frac{1}{R} \sum_{r=1}^R\left({\hat{\pi_{c}}}^{(r)}-\pi_{c }\right)^2}, \end{align*} $$

where R represent the replication times, and $\hat {.}^{(r)}$ represent the estimated value of a parameter in the rth replication. In addition, we provide the 95% highest probability density interval (HPDI) of $\boldsymbol {s}$ , $\boldsymbol {g}$ , and $\boldsymbol {\pi }$ , which is the narrowest interval that contains 95% of the posterior probability, to access the uncertainty related to these parameter estimation.

To evaluate the recovery of attribute profile estimation, we adopted the pattern-wise agreement rate (PAR), $\text {PAR}=\frac {1}{N}\sum _{i=1}^{N}\mathcal {I}(\hat {\boldsymbol {\alpha }}_i=\boldsymbol {\alpha }_i)$ , and the attribute-wise agreement rate (AAR), defined as $\text {AAR}(k)=\frac {1}{N}\sum _{i=1}^{N}\mathcal {I}(\hat {\alpha }_{ik}=\alpha _{ik}).$ Here, $\hat {\boldsymbol {\alpha }}_i$ is the estimated value of $\boldsymbol {\alpha }_i$ , and $\hat {\alpha }_{ik}$ is the estimated value of $\alpha _{ik}$ . Higher values of PAR and AAR indicate better recovery of attribute profile.

Third, to assess the recovery of Q-matrix, we calculate the Q-matrix recovery rate for each dimension ( $\text {QRR}$ ) and the average Q-matrix recovery rate ( $\text {AQRR}$ ) across the Q-matrix as follows:

$$ \begin{align*} \text{QRR(k)}=\frac{1}{R}\sum_{r=1}^R\left(1-\frac{\sum_{j=1}^J\left|\hat{q}_{jk}^{(r)}-q_{jk}\right|}{J}\right), \text{AQRR}=\frac{1}{K} \sum_{k=1}^K \text{QRR(K)}, \end{align*} $$

where $r=1,2, \ldots , R$ and $\hat {q}_{jk}^{(r)}$ is the estimated value of $q_{jk}$ in the rth replication. The higher value of AQRR and ORR indicates better recovery of $\text {Q}$ -matrix. Additionally, we also report the number of times the overall $\text {Q}$ -matrix was correctly estimated across 50 replication.

Finally, when evaluating the recovery of the hierarchical structure, we considered four criteria: the true positive rate (TPR) and the true negative rate (TNR), based on the adjacency matrix ( $\text {G}$ -matrix), as well as the reachability TPR (RTPR) and the reachability TNR (RTNR), based on the reachability matrix (R-matrix):

$$ \begin{align*}\operatorname{TPR}=\frac{\#\operatorname{True~positive~edge} }{\#\operatorname{True~positive~edge+False~ negative~ edge}} , \end{align*} $$

$$ \begin{align*}\operatorname{TNR}=\frac{\#\operatorname{True~negative~edge} }{\#\operatorname{True~negative~edge+False~ positive~ edge}} , \end{align*} $$

$$ \begin{align*}\operatorname{RTPR}=\frac{\#\operatorname{True~reachability~positive~edge} }{\#\operatorname{True~reachability~positive~edge+False~reachability~ negative~ edge}} , \end{align*} $$

$$ \begin{align*}\operatorname{RTNR}=\frac{\#\operatorname{True~reachability~negative~edge} }{\#\operatorname{True~reachability~negative~edge+False~reachability~ positive~ edge}}. \end{align*} $$

TPR is the proportion of positive edges that are correctly estimated, TNR denotes the proportion of negative edges that are correctly estimated, RTPR measures the proportion of reachability positive edges that are correctly estimated, and RTNR refers to the proportion of reachability negative edges that are correctly estimated. The higher and closer to 1 value of TPR, TNR, RTPR, and RTNR indicate a better recovery of the hierarchical structure $\text {G}$ . The number of times the overall $\text {G}$ was correctly estimated across 50 replications was also reported.

5.4 Simulation results

5.4.2 Computational time

We report the average computation time (in seconds) across all replications under each simulation condition in Table 3. Algorithm was operated in R programming language (R Core Team, 2017) on a desktop computer equipped with Intel(R) Core(TM) i5-10400 CPU @ 2.90 GHz, 16 GB RAM. The two R packages “Rcpp” (Eddelbuettel & François, Reference Eddelbuettel and François2011) and “RcppArmadillo” (Eddelbuettel & Sanderson, Reference Eddelbuettel and Sanderson2014) were used to enhance the computational efficiency. From the result in Table 3, when $N=500$ , it takes within 6 minutes across all conditions, and it takes within 19 minutes across all conditions when $N=1,000$ .

Table 3 Average computation time of the proposed algorithm under each condition in seconds

5.4.3 The impact of mini-batch sample size and cut-off value of $\text {G}$

We first explore how different levels of mini-batch size and cut off values of $\text {G}$ impact the performance of the proposed Gibbs algorithm. The recovery results of the $\text {Q}$ -matrix and the hierarchical structure $ \text {G} $ under different mini-batch sample sizes and cut-off values are documented in Appendix B of the Supplementary Material. Here, we summarize the optimal cut-off value and mini-batch size under various simulation conditions in Table 4. In this table, each cell lists the mini-batch size and cutoff value that resulted in the best recovery of $\text {Q}$ and $\text {G}$ in a specific simulation condition. For example, for the attribute hierarchy $\text {G}_1$ and sample size is 500 (test length is 20), the mini-batch size of 256, and the cut off value of 0.2 or 0.3 result in the highest recovery of Q and the mini-batch size of 256 and the cut off value of 0.3 resulted in the best recovery of $\text {G}_1$ . From this table, we can see that for the cut-off value, under the condition $ N = 500 $ , the highest recovery accuracy for both $ \text {Q} $ and $ \text {G} $ occurs when the cut-off value is 0.2, while for $ N = 1,000 $ , both 0.2 and 0.3 yield high accuracy. For the mini-batch sample size, when $N=500$ , the highest recovery accuracy for both $ \text {Q} $ and $ \text {G} $ occurs when $ N_1 = 256 $ . Among all eight conditions under $N=500$ , the optimal mini-batch sample size of $ N_1 = 128 $ appears only once, indicating that $ N_1 = 256 $ is clearly superior. When $ N = 1,000 $ , $N_1 = 512 $ outperforms $ N_1 = 256 $ only two times, indicating that the performance of $ N_1 = 256 $ and $ N_1 = 512 $ is similar. We speculate that this is because, with larger sample sizes, even if the mini-batch sample size is not very large, the continual sampling process allows the mini-batch samples to cover all the sample information since the ergodicity, leading to more accurate estimates. Based on these findings, we use mini-batch size $ N_1 = 256 $ when $ N = 500 $ and $ N_1 = 512 $ for $ N = 1,000 $ , and cut off value of 0.2 in the subsequent analysis.

Table 4 Summary of the optimal mini-bath size and cut-off value

5.4.4 Recovery of ${G}$

Table 5 displays the recovery of the hierarchical structure $\text {G}$ . Based on the results, we find that when $N=500 $ , the proportion of accurately estimating $\text {G}$ exceeds 80% across all replications and conditions, and this proportion increases to over 88% when $N=1,000$ . The TPR and TFR results indicate that, across all simulation conditions, the proposed algorithm achieves over 90% accuracy, except under $\text {G}_1$ in unbalanced condition, in estimating the hierarchical structure edges between attributes, and the RTPR and RTFR results show that the accuracy in estimating direct or indirect dependencies between attributes is even higher, reaching over 95%. Moreover, Table 6 presents the average posterior sample mean of $\text {G}$ ( $\hat {\text {G}}$ ) across 50 replications under the condition $N = 1,000$ , $J = 20$ , where the bold numbers indicate the true edges. The results show that our algorithm yields estimated values above 0.4 for the true edges between attributes. For non-existing edges but with indirect relationships between attributes, the estimated values range between 0.2 and 0.4. For edges that do not exist and have no indirect relationships between attributes, the estimated values are close to 0. This finding indicates that our method performs remarkably well in identifying the reachability relationships among attributes.

Table 5 Evaluate the recovery accuracy of the adjacency matrix $\text {G}$

Note: TPR and TFR are the true positive rate and true negative rate based on G-matrix; RTPR and RTFR are the true positive rate and true negative rate based on reachability matrix (R-matrix). The last column reports the number of times the entire G-matrix was recovered correctly.

Table 6 Estimated adjacency matrices $\hat {\text {G}}$ of the proposed algorithm for $N=1,000$ , $J=20$ condition

Note: The bolded entries in the matrix represent the real edges.

In cases where the hierarchical structure was incorrectly estimated, we found that the primary sources of error were the misestimation of directional relationships among attributes and, in some instances, the failure to identify certain edges. Notably, the recovery of structure $\text {G}_1$ was less ideal compared to the other structures when the sample size was 500. Upon further examination, we observed that most of the estimation failures were caused by incorrect directionality within the attribute hierarchy. For instance, in one replication, the true structure $1 \rightarrow 2 \rightarrow 3 \rightarrow 4$ was incorrectly estimated as $1 \rightarrow 2 \rightarrow 4 \rightarrow 3$ , in which the direction between attributes 3 and 4 was reversed. Similar issues were found under structure $\text {G}_3$ in the balanced condition, where the true edges $1 \rightarrow 2$ , $1 \rightarrow 3$ , and $1 \rightarrow 4$ were sometimes incorrectly estimated as $2/3/4 \rightarrow 1$ . Moreover, under structure $\text {G}_3$ , many incorrect estimations involved missing edges connected to attribute 1. Regardless of which edge was missing, this consistently led to the lowest QRR value for the first attribute.

In addition, the overall recovery accuracy was the lowest under the linear structure $\text {G}_1$ . To further investigate this, we conducted an additional analysis of the simulation data and found two main contributing factors: (1) a small proportion of items measuring certain attributes in the $\text {Q}$ -matrix and (2) an unbalanced distribution of examinees across the permissible attribute mastery patterns. Specifically, under both balanced and unbalanced conditions, estimation failures often occurred when fewer than 15% of items measured the fourth attribute. Moreover, in the unbalanced condition with $N = 500$ , $J=20$ , the examinees with mastery patterns $(1,0,0,0)$ and $(1,1,0,0)$ accounted for less than 10%, sometimes even below 5%, leading to insufficient information for accurate estimation. When the sample size increased to $N = 1,000$ , this issue was largely mitigated. These findings indicate that both limited item coverage for certain attributes and an unbalanced distribution of attribute mastery patterns can substantially reduce estimation accuracy.

Finally, an interesting finding is that, although the recovery of $\text {G}$ shows similar overall accuracy under both balanced and unbalanced attribute profile conditions, there are subtle differences. For example, under a sample size of 500 and test length of 20, the recovery rate under the balanced condition is higher for all structures except $\text {G}_3$ . However, when the sample size increases to 1,000 and test length increases to 40, the recovery rate under the unbalanced condition becomes higher. Additionally, under the balanced attribute profile condition, increasing the sample size and test length does not necessarily improve the recovery of $\text {G}$ . In contrast, under the unbalanced attribute profile condition, the increase in sample size and test length lead to better recovery of $\text {G}$ .

We conjecture that under the balanced condition, the sample of each attribute profile is already relatively well-balanced even at smaller sample sizes, leading to relatively good recovery rates. Therefore, increasing the sample size from 500 to 1,000 does not result in a significant improvement. In the unbalanced condition, however, the imbalance in the sample of each attribute profile results in relatively poor estimation performance at smaller sample sizes. As the sample size increases, the sample of each attribute profile increases, leading to a noticeable improvement in recovery rates.

5.4.5 Recovery of Q matrix

Table 7 presents the recovery accuracy of the $\text {Q}$ -matrix. First, across all conditions, the percentage of replications that can fully recover the whole true Q matrix, ranging from $60\%(30/50)$ to $100\% (50/50)$ , and with mean of $92\%$ . Second, overall, the proposed algorithm generally resulted in the same or higher recovery of $\text {Q}$ , calculated as the number of times the estimated $\text {Q}$ -matrix matches the true $\text {Q}$ -matrix. This improvement was observed with increasing test length and sample size across almost all conditions, except in the balanced condition with the hierarchical structure $\text {G}_2$ . In this case, the number of correct $\text {Q}$ estimations decreased from 48 to 45 when moving from $N=500$ , $J=20$ to $N=1,000$ , $J=40$ .

Table 7 Evaluate the recovery accuracy of $\text {Q}$ -matrix

Note: QRR1–QRR4 represent the correct recovery rates for the four attributes, AQRR denotes the element-wise recovery rate of the Q-matrix, and the last column reports the number of times the entire Q-matrix was recovered correctly.

From the element-wise perspective, based on the AQRR value, the element-wise recovery rate was above 97% across all simulation conditions. In the balanced conditions, under structures $\text {G}_1$ , $\text {G}_2$ , and $\text {G}_3$ , the AQRR value showed a slight decrease from $N=500$ , $J=20$ to $N=1,000$ , $J=40$ , while in other conditions, it showed a slight increase. The AQRR values are similar for sample sizes of 500 and 1,000. Similar to the overall recovery of $\text {Q}$ , the element-wise recovery rate does not always improve with an increasing sample size. We speculate that this is because, although the sample size increased from 500 to 1,000, the test length also increased from 20 to 40, meaning the number of elements in $\text {Q}$ that required estimation doubled thus increase the difficulty in estimation as well.

In addition, the recovery of $\text {Q}$ , whether considering element-wise accuracy or overall accuracy, generally shows similar results and trends under both balanced and unbalanced conditions, with only slight differences. When $N=1,000$ , the overall number of correct recoveries of $\text {Q}$ in the unbalanced condition is slightly higher than in the balanced condition.

Finally, given a sample size, the recovery of $\text {Q}$ is similar across the four considered attribute hierarchical structures. However, in terms of the total number of correct recoveries for $\text {Q}$ , the unbalanced condition with a sample size of 500 shows only 30 correct recoveries for $\text {G}_1$ (linear structure), which is noticeably lower than in the other conditions. Nevertheless, the corresponding QARR values are similar to those for the other structures. This is mainly closely related to the recovery of $\text {G}$ under that condition. Specifically, by comparing Tables 4 and 5, we found that the total number of correct estimations for $\text {Q}$ generally aligns closely with that of the hierarchical structure $\text {G}$ , with the correct estimations for $\text {Q}$ being less than or equal to those for $\text {G}$ . This suggests that if $\text {G}$ is estimated correctly, there is a high likelihood that $\text {Q}$ will also be accurately estimated. Conversely, if $\text {G}$ is estimated incorrectly, $\text {Q}$ is unlikely to be estimated correctly, since $\text {G}$ directly influences the permissible patterns for $\text {Q}$ . For unbalanced condition with a sample size of 500 under the hierarchical structure $\text {G}_1$ , only 40 out of 50 replications generate correct estimations of the overall $\text {G}$ , which is lower than those from other conditions, thus also leads to the undesired recovery of $\text {Q}$ .

5.4.6 Recovery of item parameters and attribute profile

Table 8 presents the ARMSE and ABias of the item parameters $\boldsymbol {s}$ and $\boldsymbol {g}$ , as well as the class membership parameter $\boldsymbol {\pi }$ . Since the true values of $\boldsymbol {s}$ and $\boldsymbol {g}$ are identical across all items, we further report the average HPDI (AHPDI) for $\boldsymbol {s}$ and $\boldsymbol {g}$ , defined as the mean of the lower and upper bounds of the HPDIs across all items and over 50 replications. And we also report the AHPDI for $\boldsymbol {\pi }$ over 50 replications in the figures in Appendix C of the Supplementary Material. The results indicate that the parameters $\boldsymbol {s}$ , $\boldsymbol {g}$ , and $\boldsymbol {\pi }$ exhibit good performance in terms of bias, with bias values close to zero under all simulation conditions. Across all hierarchical structures, as the sample size and test length increase from $N = 500, J = 20$ to $N = 1,000, J = 40$ , the average RMSE of $\boldsymbol {s}$ and $\boldsymbol {g}$ shows a declining trend, indicating that the accuracy of item parameter estimation improves with larger sample sizes. Simultaneously, the 95% AHPDI for parameters $\boldsymbol {s}$ , $\boldsymbol {g}$ , and $\boldsymbol {\pi }$ becomes narrower, further demonstrating that parameter estimation accuracy benefits from larger samples. Furthermore, Table 9 reports the recovery accuracy of the attribute profiles $\boldsymbol {\alpha }$ . The values of AAR and PAR also increase as the sample size grows, suggesting that the recovery accuracy of $\boldsymbol {\alpha }$ improves with larger sample sizes.

Table 8 Evaluate the estimation accuracy for the item parameters $\boldsymbol {s}$ , $\boldsymbol {g}$ , and the class membership probability parameter $\boldsymbol {\pi }$

Note: ARMSE, ABias, and AHPDI represent the average RMSE, average bias, and average high density posterior interval across all items or all latent classes.

Table 9 Evaluate the estimation accuracy for the attribute profile $\boldsymbol {\alpha }$

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

6.1 Analysis procedures

Given the total sample size, we chose three mini-batch sample sizes $N_1=512$ , 1,024, and 1,500 to implement the proposed algorithm. A total of four Markov chains were conducted for each selected batch size, and the chain length was set as 20,000 and burn in value is 10,000. In addition, the prior for each parameter is consistent with the settings used in the simulation study. Our purpose is to investigate whether the proposed algorithm can successfully recover the linear hierarchical structure $3\rightarrow 2\rightarrow 1$ among the three skills in the ECPE test, as identified in previous research (Templin & Bradshaw, Reference Templin and Bradshaw2014; Wang & Lu, Reference Wang and Lu2021).