Bayesian Estimation of Attribute Hierarchy for Cognitive Diagnosis Models

Abstract

Attribute hierarchy, the underlying prerequisite relationship among attributes, plays an important role in applying cognitive diagnosis models (CDM) for designing efficient cognitive diagnostic assessments. However, there are limited statistical tools to directly estimate attribute hierarchy from response data. In this study, we proposed a Bayesian formulation for attribute hierarchy within CDM framework and developed an efficient Metropolis within Gibbs algorithm to estimate the underlying hierarchy along with the specified CDM parameters. Our proposed estimation method is flexible and can be adapted to a general class of CDMs. We demonstrated our proposed method via a simulation study, and the results from which show that the proposed method can fully recover or estimate at least a subgraph of the underlying structure across various conditions under a specified CDM model. The real data application indicates the potential of learning attribute structure from data using our algorithm and validating the existing attribute hierarchy specified by content experts.

Keywords

cognitive diagnosis models directed acyclic graphs hierarchical structure

1. Introduction

Cognitive diagnosis models (CDMs) have been widely applied to the field of educational assessment, psychiatric diagnosis, and other social sciences (Chen et al., 2015; de la Torre & Douglas, 2004, 2008; Sorrel et al., 2016; K. K. Tatsuoka, 1984; C. Tatsuoka, 2002; von Davier, 2008). In conjunction with diagnostic assessments, this family of structured latent class models uses subjects’ observed responses to specifically designed diagnostic items to determine fine-grained classification of the underlying latent attribute patterns. However, an unsolved problem for applying CDM in a learning context is that there is still a lack of efficient statistical tools to directly learn attribute hierarchy from the observed data. An attribute hierarchy in a given domain summarizes the relationships among skills, indicating the mastery of an attribute is a prerequisite to the mastery of another one. Different attribute hierarchies have been observed in many learning contexts, such as mathematics learning (Leighton & Gierl, 2007), critical reading (Wang & Gierl, 2011), and syllogistic reasoning (Leighton et al., 2004). Research has indicated that an attribute hierarchy has a potential to play an important role in designing an effective diagnostic assessment and an adaptive learning system (Chen et al., 2018).

Most studies of attribute hierarchy within CDM assume attribute hierarchy is known and develop a series of models that incorporates the attribute hierarchy (Leighton et al., 2004; K. K. Tatsuoka, 2009; Templin & Bradshaw, 2014) or investigate the impact of attribute hierarchies on model fit and item parameter estimation (de la Torre et al., 2010) and attribute pattern classification accuracy (Tu et al., 2019). Results from these studies demonstrate that if the attribute hierarchy is correctly specified, then with the restriction on the permissible attribute patterns given a specified structure, the number of latent attribute patterns can be reduced from the full space and this can improve the model fit and item parameter estimation. Furthermore, the sample size needed to obtain stable parameter estimates from CDM calibrations will decrease as well, and this can result in faster computing time especially when there are a large number of attributes and items (Tu et al., 2019).

However, identification and validation of attribute hierarchy are usually difficult. The attribute hierarchy of a given learning context is typically determined by the domain experts based on cognitive theories. This might cause a structure misspecification problem due to subjective understandings and opinions. Few attempts in CDM literature explored the possibility of identifying attribute hierarchy from subjects’ response data. For example, several studies have developed hypothesis testing methods when some attribute structures have already been identified (e.g., Ma & Xu, 2021; Templin & Bradshaw, 2014). However, these approaches cannot be applied when no prior information of attribute hierarchy is available. Wang and Lu (2021) developed regularization-based approaches to infer attribute hierarchy directly from data. Similarly, Ma et al. (2022) proposed a penalized likelihood approach to estimate the Q matrix and restricted attribute profiles. The structure is then inferred by partial orders. These penalized likelihood approaches infer hierarchy structures through restricted attribute profile spaces, but the estimation of the actual hierarchy structure requires stepwise reconstruction of the binary representations of the partial orders.

In viewing of the aforementioned limitations of current approaches for estimating attribute hierarchy, there is considerable scope of future research on developing new methods in learning attribute hierarchy from data. Borrowing the idea of structure learning problem in the field of machine learning, we view the underlying attribute hierarchy as a structure of directed acyclic graph (DAG) and then aim to develop new Bayesian estimation approaches for learning attribute hierarchy from data within cognitive diagnosis. The rationale of framing the problem in this way is because the DAG specifies the paths connecting nodes and the distribution satisfies certain conditional independence implied by the Markov assumption. These properties are similar to the assumptions for attribute hierarchy structures in cognitive diagnosis. However, different from traditional structure learning problems, which consider the structures of observational data without hidden variables, or the hidden variables are not included in the structures, CDMs are restricted latent class models with constraints on model parameters, and the structure of latent variable is the vital information to be learned for CDM with attribute hierarchies. Therefore, the existing algorithms of structure learning (e.g., Colombo & Maathuis, 2014; Colombo et al., 2012) cannot be directly applied to our setup. Thus, the goal of this study is to develop an efficient algorithm to estimate the hierarchy structure along with the estimation of CDM model parameters and subjects’ latent attribute patterns. Our proposed algorithm explores the sample space of DAGs through random walks on transitive reductions and is applicable to a general class of CDMs under the Bayesian framework. In addition, the proposed algorithm gives exploratory discovery of hierarchical structures with posterior probabilities of potential DAGs (or edges in DAG); the posteriors provide more information on the possible prerequisite relationships among attributes than a strict estimated structure does.

The remainder of this article is organized as follows. Section 2 introduces the general CDM structures and proposes the Bayesian framework for attribute hierarchy. Section 3 introduces the developed Metropolis within Gibbs algorithm for structure learning and model estimation. Section 4 reports the simulation results of the proposed sampling algorithm using a simple CDM model as an example. Section 5 applies the algorithm to a real-world data. Finally, Section 6 provides a summary of our method and findings, discusses the implications of the study, and provides potential future research directions.

2. Bayesian Formulation for CDMs With Attribute Hierarchy

In this section, we first provide a general setup for CDM and then introduce the Bayesian estimation framework for CDMs with attribute hierarchy. Consider a scenario that N subjects complete J diagnostic questions measuring K attributes. The binary response vector, indicating the response to each item is right or wrong, from subject i is denoted by $Y_{i} = (Y_{i 1}, \dots, Y_{i J})$ . Within the CDM framework, a latent attribute profile with K binary components is assumed for each subject, generically denoted by α. That is, $α \in {0, 1}^{K}$ , indicating the nonmastery or mastery of K attributes in an educational learning context. Under the independent structure, the $2^{K}$ attribute profile patterns are all permissible, and they are represented by the index set $C = {0, 1, 2, \dots {,2}^{K} - 1}$ . Additionally, CDMs use a $J \times K$ binary $Q = {q_{j}}_{j = 1}^{J}$ matrix to indicate whether an attribute is measured by an item (K. K. Tatsuoka, 1985). Here, q _j is the jth row of Q, implying the required attributes for item j. CDMs then impose assumptions about how each subject’s latent attribute profile and each item’s characteristics together determine the subject’s correct response probability to that item. Specifically, in the general setup of CDMs (e.g., de la Torre, 2011; von Davier, 2008), $Y_{i j}$ s are modeled as independent Bernoulli variables with parameter $θ_{i j} = f (μ_{i j})$ , that is

p (Y_{i j} = 1 | α_{i}, q_{j}, β_{j}) = f (μ_{i j}) .

Here, f is a link function, and $μ_{i j} = g (β_{j}, α_{i})$ for item j with parameter $β_{j}$ and subject i with attribute profile $α_{i}$ . $μ_{i j}$ can be formulated as

\begin{array}{l} μ_{i j} = g (β_{j}, α_{i}) = β_{j 0} + β_{j 1} α_{i 1} + β_{2} α_{i 2} + \dots + β_{j K} α_{i k} + \sum_{k < k^{'}} β_{j k k^{'}} α_{i k} α_{i k^{'}} \\ + \sum_{k < k^{'} < k^{*}} β_{j k, k^{'}, k^{*}} α_{i k} α_{i k^{'}} α_{i k^{*}} + \dots + β_{j 1, 2, \dots, K} \prod_{k = 1}^{K} α_{i k} . \end{array}

Note that the parameter vector indexed by $β_{j} = (β_{j 0}, β_{j 1}, \dots, β_{j k}, \dots, β_{j k, \dots, k^{'}}, \dots, β_{j 1, \dots, K}), k < \dots < k^{'}$ is equivalent to $β_{j} = {β_{j c}}_{c = 0}^{2^{K} - 1}$ . The choice of the link function and the sparsity (zeros and nonzeros) of $β_{j}$ are specified by the assumption of a given CDM and q _j (e.g., Chen et al., 2020; de la Torre, 2011).

The existence of an attribute hierarchy directly impacts the cardinality of the attribute profile space, thus influencing $p (Y_{i j} = 1 | α_{i}, q_{j}, β_{j})$ . To start with, we first formulate an attribute hierarchical structure as a DAG. This is because the prerequisite relationship is single-directed, and two attributes cannot be each other’s prerequisite at the same time. Figure 1 shows four types of hierarchy structure proposed by Leighton et al. (2004). Each of these structures is clearly a DAG with attributes as nodes and a directed edge from k to $k^{'}$ represents that $α_{k}$ is the prerequisite of $α_{k^{'}}$ . From now on, we use $G (V, E)$ to denote an attribute hierarchy, where V is the set of nodes (latent attributes) and E is the set of directed edges (prerequisite relationships), and we use the notation G to represent the adjacency matrix of structure $G (V, E)$ and thus simplify the hierarchical structure as G for the remainder of this article.

Figure 1.

Illustration of four types of attribute hierarchy structure. (a) Linear. (b) Convergent. (c) Divergent. (d) Unstructured.

The attribute profile space is determined by a specific G. For example, for the linear structure presented in Figure 1a, an attribute profile with $α_{2} = 1$ and $α_{1} = 0$ is not permissible as this violates the assumption that attribute 1 is the prerequisite for attribute 2. Thus, the cardinality of the attribute profile space given a G will be smaller than $2^{K}$ , which is the largest value when we do not assume any specific hierarchical structure among attributes. Another component one needs to carefully specify is the Q matrix when a G is assumed. In fact, there is a debate in the literature regarding whether Q matrix needs to reflect the specified attribute hierarchy. On the one hand, several research studies assume items do not necessarily reflect the specified structure, thus have used an unstructured item-by-attribute Q matrix in their studies regarding attribute hierarchy in CDM (e.g., de la Torre et al., 2010; Templin & Bradshaw, 2014). The unstructured Q matrix only labels the attribute(s) directly measured by an item as 1(s). On the other hand, another group of researchers think a structured Q matrix that reflects the specified attribute hierarchy is necessary (e.g., Leighton et al., 2004; Tu et al., 2019). The structured Q matrix not only labels the attributes directly measured by an item but also their prerequisites.

In our study, we assume the Q matrix is known, and the goal is to estimate the latent attribute profiles, item parameters, and the DAG directly from the observed response data. Thus, we prefer to assume the Q matrix is unstructured, as the true attribute hierarchy is unknown. In addition, another benefit of using an unstructured Q is that the sparsity of model parameters $β$ is fully determined by Q and a CDM model assumption. In other words, given a specified CDM assumption (i.e., conjunctive, disjunctive, or compensatory assumption), the sparsity of $β$ has a one-to-one relationship with the unstructured Q (Chen et al., 2020), and this helps simplify our model estimation procedure. In addition, even though the Q matrix does not reflect the specified attribute hierarchy structure, the estimation method described in the following section already restricts the possible structure through the latent attribute profile space, so that the estimation result from our method will still align with the implied structure if there were any.

In summary, let $B = (β_{1}, \dots, β_{J})$ be the collection of item parameters for J items. With the conditional independence assumption, the likelihood of the response vector of N subjects, $Y = (Y_{1},..., Y_{N})^{Τ}$ , is

p (Y | A, Q, B, G) = \prod_{i} \prod_{j} p (Y_{i j} | α_{i}, q_{j}, β_{j}, G),

where $A = (α_{1}, \dots, α_{N})^{Τ}$ is a $N \times K$ matrix representing attribute profiles of all subjects. As noted above, due to the fact that the information about the Q-matrix is implied by the sparsity of B, the likelihood (2) can be simplified as $p (Y | A, B, G) = \prod_{i} \prod_{j} p (Y_{i j} | α_{i}, β_{j}, G)$ . To establish the Bayesian estimation framework for G, let $p (G)$ be the prior for all possible DAGs. The posterior distribution of G is

p (G | Y) = \frac{\int_{B} (\sum_{A} p (Y | A, B, G) p (A, B, G)) d B}{p (Y)} .

The derivation of the posterior of G given Y is nontrivial, and we use the Bayesian directed graphical models for factorization of the joint probabilities of p(A, B, G) and propose the following two tasks to sequentially update G, $A$ , and $B$ .

3. Metropolis Within Gibbs Sampling Algorithm

In this section, we first focus on the estimation of $A$ , $B$ , and G for a general class of CDMs formulated by Equation 1. We will then illustrate the proposed algorithm using a specific CDM as an example. To start with, the CDMs with attribute hierarchy can be formulated by a graphical model as shown in Figure 2.

Figure 2.

Graphical model of the Bayesian framework for cognitive diagnosis model with hierarchy structure.

3.1. Gibbs Sampling for Model Parameters Under Hierarchy Structure

We assume that the hierarchy structure G only affects the attribute profiles proportions $π$ and hence the profiles A. Then, the problem of attribute hierarchy estimation under this Bayesian graphical model framework essentially depends on $p (G | Y, A, B, π)$ . The full conditional distribution of G is as follows:

p (G | Y, A, B, π) \propto p (Y | A, B) p (A | π, G) p (π | G) p (G) .

According to (4), we propose to sequentially update A, $π$ , B, and G based on their posteriors through the Gibbs sampling algorithm. Specifically, we assume $π$ follows a Dirichlet distribution $D i r (δ_{0})$ , where $δ_{0}$ is a $2^{K}$ hyperparameter vector, and A follows a multinomial distribution with parameters $π$ .

Algorithm 1 shows the Gibbs sampling algorithm for all model parameters and the hierarchy structure illustrated in Figure 2. Algorithm 1 mainly consists of the procedures that update CDM parameters given G (Steps 3–5) and that update G given a CDM model (Step 6). We illustrate these procedures in the following sections.

Algorithm 1 Gibbs Sampling Algorithm for Attribute Hierarchy Estimation.
1: Initialize with an non-hierarchical structure $G^{(0)} = 0$ , attribute pattern proportions $π_{0}$ and individual attribute profiles A ⁰, and model parameters $B^{(0)}$ . 2: for t in $1 : T$ do 3: Update $A^{(t)}$ from $p (A \| Y, G^{(t - 1)}, π^{(t - 1)}, B^{(t - 1)})$ $\propto p (Y \| A, B^{(t - 1)}) p (A \| π^{(t - 1)}, G^{(t - 1)}) p (π^{(t - 1)} \| G^{(t - 1)})$ . 4: Update $π^{(t)}$ from $p (π \| A^{(t)}, G^{(t - 1)}) \propto p (A^{(t)} \| π, G^{(t - 1)}) p (π \| G^{(t - 1)})$ . 5: Update $B^{(t)}$ from $p (B \| Y, A^{(t)}) \propto p (Y \| B, A^{(t)}) p (B)$ . 6: Update $G^{(t)}$ from $p (G \| Y, A^{(t)}, π^{(t)}, B^{(t)}) \propto p (Y \| B^{(t)}, A^{(t)}, G) p (A \| π^{(t)}, G) p (π^{(t)} \| G) p (G)$ . 7: end for

Algorithm 1 Gibbs Sampling Algorithm for Attribute Hierarchy Estimation.

1: Initialize with an non-hierarchical structure

G^{(0)} = 0

, attribute pattern proportions

π_{0}

and individual attribute profiles A ⁰, and model parameters

B^{(0)}

.
2: for t in

1 : T

do
3: Update

A^{(t)}

from

p (A | Y, G^{(t - 1)}, π^{(t - 1)}, B^{(t - 1)})

\propto p (Y | A, B^{(t - 1)}) p (A | π^{(t - 1)}, G^{(t - 1)}) p (π^{(t - 1)} | G^{(t - 1)})

.
4: Update

π^{(t)}

from

p (π | A^{(t)}, G^{(t - 1)}) \propto p (A^{(t)} | π, G^{(t - 1)}) p (π | G^{(t - 1)})

.
5: Update

B^{(t)}

from

p (B | Y, A^{(t)}) \propto p (Y | B, A^{(t)}) p (B)

.
6: Update

G^{(t)}

from

p (G | Y, A^{(t)}, π^{(t)}, B^{(t)}) \propto p (Y | B^{(t)}, A^{(t)}, G) p (A | π^{(t)}, G) p (π^{(t)} | G) p (G)

.
7: end for

3.2. Efficiently Sampling of G

The main challenge of Algorithm 1 is in developing a sampling scheme to generate $G s$ from $p (G | Y, A^{(t)}, π^{(t)}, B^{(t)})$ , which may not be in a closed form. With a constrain of a specific structure and a reasonable value of K, the cardinality of G is generally not extremely large, and we proposed to use a Metropolis–Hastings (M–H) algorithm to sample G from the whole space of all possible DAGs through adding or removing one directed edge at each iteration. A key step in this M–H algorithm is to propose a transition step to ensure an efficient search of the candidate sample space. We first introduce the proposed transition step and then present the developed M–H algorithm.

An efficient transition step

To start with, DAGs with different numbers of edges might have the same dependent relationships. This is known as Markov equivalence property of DAGs (Chickering, 2002). For example, Figure 3 shows two DAGs that have the same directed connection (the so-called reachability). The extra edge in Figure 3b shows the connection of $1 \to 4$ , which is already implied by the path of $1 \to 2 \to 4$ . In terms of conditional independence, under the global Markov condition, both graphs imply that Node 3 is independent of Nodes 2 and 4 given Node 1.Therefore, Figure 3a with a simpler structure can represent both DAGs, and a random walk from Fig. 3(a) to (b) is inefficient. To make efficient transitions, we propose a graph generating mechanism using the notion of transitive reduction to construct a simple structure. A transitive reduction is a subgraph of G with the fewest edges that maintain the same reachability as G, and it can be retrieved from any given DAGs by selecting edges that are the only paths between its endpoints (Aho et al., 1972). In this regard, the subgraph is unique for DAGs in the same equivalence class. Adopting this idea, we propose to only allow transitions between transitive reduction DAGs.

Figure 3.

Illustration of directed acyclic graphs belonging to the same Markov equivalence class.

Proposed M–H Sampler

Based on the concept of transitive reductions, we developed an M–H method, documented by Algorithm 2, to sample $G s$ from $p (G | Y, A^{(t)}, π^{(t)}, B^{(t)})$ (line 6 of Algorithm 1). Specifically, starting with any transitive reduction G, with a prespecified probability p, a new directed edge is added between two nodes that are not connected through any paths. This guarantees the newly added edge is the only path between its endpoints; thus, the newly proposed $N = 500$ remains a transitive reduction. Otherwise, an existing edge is removed. Note that if the current $G^{*}$ is a transitive reduction, then the newly proposed $G^{*}$ after removing one edge is still a transitive reduction, because the remaining edges are still the only paths between its endpoints.

Let $T (G, G^{*})$ denote the transition probability from G to the proposed new $G^{*}$ , then

T (G, G^{*}) = \frac{1}{2 \times number of pairs of non - connected nodes},

if the proposed step is adding an edge,

T (G, G^{*}) = \frac{1}{number of existing edges}, if the proposed step is removing an edge .

Correspondingly the acceptance ratio of a proposed M–H sampler $G^{*}$ is

r = min {1, \frac{f (G^{*}) T (G^{*}, G)}{f (G) T (G, G^{*})}},

= min {1, \frac{p (Y | G^{*}, A, B) p (A | π, G^{*}) p (π | G^{*}) p (G^{*}) T (G^{*}, G)}{p (Y | G, A, B) p (A | π, G) p (π | G) p (G) T (G, G^{*})}},

where $p (G)$ is assumed to be uniform among possible equivalent classes of DAGs, and $π | G^{*}$ follows a Dirichlet distribution. When G is a graph with no edges, which implies an independent structure, $π$ follows the original prior of Dirichlet distribution with parameters $δ_{0}$ of length $2^{K}$ . When G imposes a specific hierarchical structure, we use the aggregation property of Dirichlet distribution, such that $π$ now follows a Dirichlet distribution with parameters $δ_{G}$ of length $| C_{G} |$ , where C_G is the index set of permissible α given G, $| C_{G} | < | C |$ . In particular, we have

p (π | G) \propto \prod_{c = 1}^{| C_{G} |} {(\sum_{\tilde{c}} π_{\tilde{c}})}^{\sum_{\tilde{c}} δ_{\tilde{c}}},

where $\tilde{c} \in C$ represents the index of attribute profile patterns that are coerced into the permissible pattern c under the structure G. The computation of coerced permissible patterns is discussed in Section 3.3. For computation efficiency, we collapsed A in the M–H samplers, such that

p (Y | G, π, B) = \prod_{i = 1}^{N} \sum_{c = 1}^{| C_{G} |} p (Y_{i} | G, α_{i} = α_{c}, B) π_{c} .

The acceptance ratio then becomes

r = min {1, \frac{p (Y | G^{*}, π, B) p (π | G^{*}) T (G^{*}, G)}{p (Y | G, π, B) p (π | G) T (G, G^{*})}} .

The proposed transition step in Algorithm 2 is irreducible, because from any transitive reduction DAG G in the Markov equivalence class, one can reach a new transitive reduction DAG $G^{*}$ by removing the edges in G sequentially and adding the edges in $G^{*}$ one by one thereafter. This irreducibility of the proposed transition step guarantees Algorithm 2 can search through the whole space of Markov equivalent DAGs. Another remark is on the prespecified probability p to control the step of proposing different structures. We hope the algorithm to explore the sampling spaces of other model parameters conditioned on a given structure and thus propose to change the structure after a few iterations of exploration. The prespecified p can be set to 1, where at each iteration, the algorithm proposes a new structure. We suggest to select p based on the acceptance ratio of the proposed samples in the Metropolis step, where the acceptance ratio between 20% and 30% is preferred.

Algorithm 2 Metropolis–Hasting Algorithm for Sampling Attribute Hierarchy Structure.
1: The current state of the structure is G and attributes profiles $A$ , with a pre-specified probability p 2: General a u from Uniform $(0, 1)$ : 3: if $u < p$ then 4: Propose a new $G^{}$ by randomly adding one directed edge between two nodes that are not connected by any paths. 5: else if $u > p$ then 6: Propose a new $G^{}$ by randomly removing one existing edge. 7: end if 8: Compute the acceptance ratio in Equation (9) and generate w from Uniform $(0, 1)$ : 9: if $w < r$ then 10: Accept the new $G^{*}$ . 11: else 12: Stay with the current G. 13: end if

Algorithm 2 Metropolis–Hasting Algorithm for Sampling Attribute Hierarchy Structure.

1: The current state of the structure is G and attributes profiles

A

, with a pre-specified probability p
2: General a u from Uniform

(0, 1)

:
3: if

u < p

then
4: Propose a new

G^{*}

by randomly adding one directed edge between two nodes that are not connected by any paths.
5: else if

u > p

then
6: Propose a new

G^{*}

by randomly removing one existing edge.
7: end if
8: Compute the acceptance ratio in Equation (9) and generate w from Uniform

(0, 1)

:
9: if

w < r

then
10: Accept the new

G^{*}

.
11: else
12: Stay with the current G.
13: end if

3.3. Sampling CDM Model Parameters

The derivation of full conditional distributions for A, $π$ , and B is explicit for a general class of CDMs following the Bayesian formulation (e.g., Chen et al., 2020). For example, for DINA model, these posteriors can be derived explicitly using certain specified priors, such as the formulations in Culpepper (2015). Even if the derivation of posteriors is explicit, the computation of $p (A | π)$ and $p (π | G)$ in each iteration needs to be sequentially adapted. This is because different G structures imply different constrained spaces of A and $π$ . The A at the current iteration may not comply with the space implied by a new $G^{*}$ in the next one.

When an attribute pattern $α_{c}$ under independent structure, indexed by $c \in C$ , is not permissible under the proposed structure G, its corresponding permissible pattern should comply with the constraint G. Denote the new pattern under structure G by $c_{G} \in C_{G}$ . By the aggregation property of Dirichlet distribution, the $π | G$ becomes a vector of length $| C_{G} |$ following an aggregated Dirichlet distribution, with

π_{c_{G}} = \sum_{α_{c} \overset{G}{=} α_{c_{G}}; c \in C} π_{c},

δ_{c_{G}} = \sum_{α_{c} \overset{G}{=} α_{c_{G}}; c \in C} δ_{c},

where $α_{c} \overset{G}{=} α_{c_{G}}$ means $α_{c}$ and $α_{c_{G}}$ are equivalent under structure G. For example, consider a simple case for $K = 3$ , and the current G has one edge from $α_{1}$ to $α_{2}$ , as shown in Figure 4a. The space of A given G is ${(0, 0, 0)$ , $(1, 0, 0)$ , $(1, 1, 0)$ , $(0, 0, 1)$ , $(1, 0, 1)$ , $(1, 1, 1)}$ . Two attribute patterns, $(0, 1, 0)$ and $(0, 1, 1)$ , are not permissible with G. Because the random walk approach only modifies one edge at each proposed step, there will be two simple scenarios of adding or removing a prerequisite relationship. If removing an edge for the next step, resulting in a new $G_{1}^{*}$ (Figure 4b); then, $α_{1}$ is no longer the prerequisite of $α_{2}$ . This implies all eight attribute patterns are permissible in the space of A given $G_{1}^{*}$ in Figure 4b. Then, $p (π | G_{1}^{*})$ is computed as the Dirichlet distribution with $2^{K}$ variables. Similarly, if the next step is to add an edge, then a new prerequisite relationship is proposed, resulting in a $G_{2}^{*}$ (Figure 4c). In this case, the space of $A | G_{2}^{*}$ is ${(0, 0, 0)$ , $(1, 0, 0)$ , $(1, 1, 0)$ , $(1, 0, 1)$ , $(1, 1, 1)}$ . That is, one permissible attribute pattern under G (Figure 4a), $(0, 0, 1)$ , is no longer permissible and the subjects under which will be reclassified into $(0, 0, 0)$ under the new $G_{2}^{*}$ in Figure 4c. This causes a reduced number of permissible patterns under $G_{2}^{*}$ . In this case, based on the aggregation property of Dirichlet distribution, we only need to collapse the impermissible attribute patterns into the permissible ones under the reduced space, and $p (π | G_{1}^{*})$ is an aggregated Dirichlet distribution with a collapsed hyperparameter.

Figure 4.

Illustration of directed acyclic graphs and permissible patterns. (a) G. (b) $G_{1}^{*}$ after removing an edge. (c) $G_{2}^{*}$ after adding an edge.

In the case that an edge is eliminated, the number of permissible mastery patterns increases and all current patterns are still permissible under the new structure; therefore, they can be directly used in the calculation for acceptance ratio (as conditioned on the current patterns). Then, in the later steps of Gibbs sampling, attribute patterns are updated through Dirichlet-multinomial distribution, where all attribute profiles are updated sequentially relying on the updated profiles only, that is, $P (α_{i} = α_{c} | G)$ is proportional to the number of subjects among the existing $i - 1$ subjects classified as $α_{c}$ . The Dirichlet parameter $π$ is then updated and conditioned on the updated attributes.

3.4. Summary of the Proposed Algorithms

In summary, incorporating Algorithm 2 into line 6 of Algorithm 1, we have the Metropolis-within-Gibbs (MH-Gibbs) algorithm for sampling the hierarchical structure and model parameters together. The algorithm is applicable to general CDMs, with explicit derivation of the posteriors of CDM model parameters. Given a specific CDM model, essentially, one only needs to modify the sampling procedures for the CDM part (Steps 3–5 in Algorithm 1) and the proposed MH algorithm for sampling G (Algorithm 2) is applicable to any CDMs. As an illustration, we present the Gibbs sampling procedures for DINA model given G here. For DINA model, $B = (s, g)$ where s is a vector of slipping parameters and g is a vector of guessing parameters. Assume the prior for $π$ is Dirichlet $(δ)$ , the priors for s and g are Beta $(a_{s}, b_{s})$ and Beta $(a_{g}, b_{g})$ , respectively. Under the proposed random walk, the prior for G is uniform on all transitive reductions of DAG of size K. Therefore, the posteriors for Gibbs samplers of A, $π$ , s, and g are

P (α_{i} = α_{c} | Y_{i}, G, s, g, α_{1 : i - 1}) \propto p (Y_{i} | G, s, g, α_{i} = α_{c}) i_{c} for c \in C_{G}, i = 1, \dots, N,

π | A, G \sim Dirichlet (δ_{G}), δ_{G} = (δ_{c} + n_{c}), n_{c} = number of subjects having α_{c}, c \in C_{G},

s_{j} | Y, A, g_{j} \sim Beta (n_{10} + a_{s}, n_{11} + b_{s}) ℐ (0 \leq g_{j} < 1 - s_{j} \leq 1),

g_{j} | Y, A, s_{j} \sim Beta (n_{01} + a_{g}, n_{00} + b_{g}) ℐ (0 \leq g_{j} < 1 - s_{j} \leq 1),

where i_c is the number of subjects among the current $_{g}$ subjects that are classified as possessing attribute profile $α_{c}$ , $n_{11}$ and $n_{10}$ are the numbers of correct and incorrect responses to item j from subjects who possess all required skills, and $n_{01}$ and $n_{00}$ are the numbers of correct and incorrect responses to item j from subjects who miss at least one required skill.

4. Simulation Study

In this section, we examined the performance of our proposed algorithm assuming response data follow the DINA model. To demonstrate the generalizability of the proposed algorithm, in addition to this simple model, we implemented the algorithm for general restricted latent class models discussed in Section 2. The results can be found in the Appendix. For this simulation study, we use different sample sizes ( $N = 500$ and 1,000), test lengths ( $J = 20$ and $J = 40$ ), numbers of attributes ( $K = 4$ and 7), attribute dependence levels ( $ρ = 0$ and $0.5$ ), and attribute hierarchy structures (shown in Figures 5 and 6). To evaluate the performance of the proposed algorithm in the case of balanced and unbalanced attribute profile distribution, we use two different distributions to simulate the true attribute profile. Specifically, when $ρ = 0$ , the profiles are randomly generated from one of its permissible pattern; thus, the true $π$ is relatively balanced among permissible patterns. When $ρ = 0.5$ , the profiles are generated from a higher order probit model that

Figure 5.

Four hierarchy structures in simulation study for K = 4. (a) G 1. (b) G 2. (c) G 3. (d) G 4.

Figure 6.

Four hierarchy structures in simulation study for K = 7. (a) G 1. (b) G 2. (c) G 3. (d) G 4.

α_{i k} = {\begin{array}{l} 1, & if Z_{i k} \geq 0 \\ 0, & otherwise \end{array},

where $Z_{i} = (Z_{i 1}, \dots, Z_{i K})^{'}$ follows a multivariate normal distribution with mean zeros and a covariance matrix $Σ = (1 - ρ) l_{K} + ρ 1 1^{'}$ , and the $2^{K}$ attribute profiles are converted into permissible patterns of the given hierarchy structure. Using this generating approach, the true $π$ is unbalanced and the proportions of all zeros or all ones are higher than others. The true slipping and guessing parameters are fixed at 0.2.

To evaluate the simulation performance, we ran 100 replications for all cases. First, we checked the convergence of the proposed algorithm using the multivariate proportional scale reduction factor (Brooks & Gelman, 1998), and the results indicate the algorithm reaches the convergence criterion after 6,000 iterations. Thus, we ran all chains of 20,000 iterations with a 15,000 burn-in period. The hyperparameter for Dirichlet prior of $π$ is set to be $δ_{0} = (1, \dots,1)$ , and the hyperparameters for slipping and guessing parameters are set to be $a_{s} = b_{s} = a_{g} = b_{g} = 1$ . All model parameters are estimated by posterior means. The hierarchy structure $\hat{G}$ is inferred by the inclusion probability on each element of the adjacency matrices from the posterior samples. The classification of attribute profiles is evaluated by overall pattern accuracy rate (PAR), which is defined by $\frac{1}{N} \sum_{i} ℐ {{\hat{α}}_{i} = α_{i}}$ and attribute accuracy rate (AAR), defined by $\frac{1}{N} \sum_{i} ℐ {{\hat{α}}_{i k} = α_{i k}}, k = 1, \dots, K$ . The estimation of model parameters is evaluated by the root mean squared error (RMSE):

{RMSE}_{s} = \sqrt{\frac{\sum_{j} {({\hat{s}}_{j} - s_{j})}^{2}}{J}}, {RMSE}_{g} = \sqrt{\frac{\sum_{j} {({\hat{g}}_{j} - g_{j})}^{2}}{J}} .

4.1. Simulation Results

Recovery of G for K = 4

For a DAG represented by an adjacency matrix, its $(i, j)$ th element represents a potential edge from attribute i to attribute j, indicating attribute i is the prerequisite of attribute j. The final estimated adjacency matrix is presented by a weighted matrix, where the $(i, j)$ th element represents the proportions of edge $i \to j$ included in the 100 replications of all posterior samples. A hierarchical structure can then be recovered from the weighted matrix with a specified threshold of inclusion probability. Table 1 reports the estimated hierarchy structures using weighted adjacency matrix under four different specified structures for $K = 4$ and $J = 20$ averaging over 100 replications. First, the estimated inclusion probabilities of all attribute hierarchy structures are comparable when $ρ = 0$ and $ρ = 0.5$ (columns 2 and 3 of Table 1), indicating the distribution of attribute profile has little impact on the estimation of G. Estimation in N = 1,000 case has a slightly higher inclusion probability for true edges than in $N = 500$ case. Second, we note that for all considered attribute hierarchical structures, the collection of estimated edges is a subset of the list of true edges. As discussed in Section 3.2, some edges in DAG can be implied from other paths, and removing those edges will not affect the hierarchical structure implied by the graph. Note that the extra edges estimated from our algorithm in Table 1 for G ₁, G ₂, and G ₄ are all implied edges, which do not violate the true structure. Therefore, based on the true structure, the proposed method correctly estimated the true edges or the implied edges. Lastly, we observe that the inclusion probabilities for true edges are around 0.5, the inclusion probabilities for implied edges are smaller than those for true edges, and the inclusion probabilities for directed edges not implied in the true structure are all 0 in $K = 4$ case. This suggests that the estimated hierarchy structures of 100 replications are either the true structure or a subgraph of the true structure. To further illustrate this, Table 2 summarizes the overall structure estimation over 100 replications for $ρ = 0$ and $ρ = 0.5$ . The threshold for inclusion probability among posterior samples is fixed as 0.5. The proposed method always identifies the true structure or a subset of it. In fact, for those estimated subset structures, the inclusion probabilities for missing edges are around 0.4, and the inclusion probabilities for nonexisting edges are 0, so one may consider a loose threshold in practice to get a better recovery of G.

Table 1.

Estimated Structures $\hat{G}$ in Adjacency Matrices of the Proposed Method for K = 4 and J = 20 Cases

G	$N = 500, J = 20, ρ = 0$				$N = 1,000, J = 20, ρ = 0$				$N = 500, J = 20, ρ = 0.5$				N = 1,000 $, J = 20, ρ = 0.5$
G1	.000	.544	.208	.153	.000	.545	.184	.183	.000	.523	.199	.159	.000	.529	.202	.141
	.000	.000	.557	.193	.000	.000	.542	.184	.000	.000	.521	.174	.000	.000	.512	.164
	.000	.000	.000	.530	.000	.000	.000	.535	.000	.000	.000	.498	.000	.000	.000	.473
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
G2	.000	.522	.544	.139	.000	.530	.529	.145	.000	.510	.516	.132	.000	.531	.522	.130
	.000	.000	.000	.497	.000	.000	.000	.499	.000	.000	.000	.490	.000	.000	.000	.485
	.000	.000	.000	.496	.000	.000	.000	.509	.000	.000	.000	.481	.000	.000	.000	.496
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
G3	.000	.497	.495	.500	.000	.538	.541	.531	.000	.485	.497	.505	.000	.534	.511	.522
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
G4	.000	.502	.542	.248	.000	.522	.545	.230	.000	.500	.502	.253	.000	.518	.533	.261
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.507	.000	.000	.000	.523	.000	.000	.000	.505	.000	.000	.000	.518
	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000	.000

Note. $\hat{G}$ is denoted by the inclusion probabilities of each entry in the sampled $K \times K$ adjacency matrices. The proportions are calculated from 100 replications of each case across all posterior samples. True edges are in bold.

Table 2.

Frequency of Corrected Estimated Structures Among 100 Replications for K = 4 and J = 20 Cases

		N = 500, J = 20			N = 1,000, J = 20
	Structure	True	Subset	False	True	Subset	False
$ρ = 0$	G1	73	27	0	77	23	0
	G2	72	28	0	72	28	0
	G3	72	28	0	76	24	0
	G4	70	30	0	72	28	0
$ρ = 0.5$	G1	54	46	0	61	39	0
	G2	65	35	0	64	36	0
	G3	56	44	0	72	28	0
	G4	60	40	0	70	30	0

Note. The threshold of inclusion probability is set to be $0.5$ . “True” refers to estimated structures (including implied edges) equivalent to the true structure, “Subset” refers to estimated structures missing some edges, and “False” refers to estimated structures with extra edges or wrong directions.

Furthermore, a better performance can be observed when both the sample size and test length increase. Specifically, N = 1,000 and $J = 40$ case gives a better performance than N = 1,000 $and J = 40$ case. Tables 3 and 4 report the estimated structures and frequency of correctly estimated structures over 100 replications when N = 1,000 and $J = 40$ . The final estimated inclusion probabilities of true edges are all greater than 0.5, and the frequency of correctly identified structures is also higher compared with the condition $N = 500 and J = 20$ .

Table 3.

Estimated Structures Ĝ_ij in Adjacency Matrices of the Proposed Method for K = 4 and J = 40 Cases

G	N = 1,000, J = 40, ρ = 0				N = 1,000, J = 40, ρ = 0.5
G1	.000	.566	.198	.145	.000	.566	.184	.142
	.000	.000	.576	.189	.000	.000	.559	.187
	.000	.000	.000	.562	.000	.000	.000	.547
	.000	.000	.000	.000	.000	.000	.000	.000
G2	.000	.535	.548	.139	.000	.521	.541	.134
	.000	.000	.000	.536	.000	.000	.000	.507
	.000	.000	.000	.534	.000	.000	.000	.528
	.000	.000	.000	.000	.000	.000	.000	.000
G3	.000	.547	.547	.532	.000	.530	.525	.535
	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.000	.000	.000	.000	.000
G4	.000	.528	.549	.249	.000	.526	.547	.222
	.000	.000	.000	.000	.000	.000	.000	.000
	.000	.000	.000	.542	.000	.000	.000	.542
	.000	.000	.000	.000	.000	.000	.000	.000

Table 4.

Frequency of Corrected Estimated Structures Among 100 Replications for K = 4 and J = 40 Cases

		N =1,000
	Structure	True	Subset	False
$ρ = 0$	G1	89	11	0
	G2	77	23	0
	G3	82	18	0
	G4	81	19	0
$ρ = 0.5$	G1	82	18	0
	G2	75	25	0
	G3	77	23	0
	G4	77	23	0

Recovery of CDM parameters for K = 4

Table 5 reports the overall attribute classification accuracy rate, the estimation accuracy rate for each attribute, and the RMSEs and biases of slipping and guessing parameters across all conditions. Again, the parameter estimation results are similar when $ρ = 0$ and $ρ = 0.5$ . The proposed method accurately classifies the attribute profiles of subjects using the estimated structure, with PAR greater than 0.957 and AAR greater than 0.982 across all conditions. We also note that for a given attribute hierarchy structure, the classification accuracy for the attribute at a lower level tends to be higher than those at a higher level. For example, for the linear stricture G ₁ that has $α_{4}$ as the lowest level and $α_{1}$ as the highest level, the AAR for $α_{4} > α_{3} > α_{2} > α_{1}$ . The RMSEs for s and g are smaller than 0.035 across all conditions, indicating an accurate estimation of item parameters for DINA. For $N = 1,000$ case, the attribute classification accuracy is slightly better than $N = 500$ case, and the RMSEs are smaller for a larger sample size. The average biases of the slipping and guessing parameters are close to 0.

Table 5.

Estimated Attribute Profile Accuracy and Item Parameter RMSEs and Biases of the Proposed Method for K = 4 and J = 20 Cases

	Structure	PAR	AAR1	AAR2	AAR3	AAR4	RMSE_s	Bias_s	RMSE_g	Bias_g
$N = 500$	G1	.978	.988	.992	0.996	1	.017	.003	.022	.006
$ρ = 0$	G2	.978	.985	.992	0.991	0.999	.017	−.008	.022	.001
	G3	.957	.986	.988	0.989	0.989	.025	.003	.020	−.001
	G4	.967	.986	.990	0.991	0.995	.021	.002	.021	.006
$N = 500$	G1	.977	.985	.994	0.997	0.999	.018	.006	.034	.001
$ρ = 0.5$	G2	.972	.983	.993	0.992	0.999	.018	.003	.033	−.003
	G3	.961	.982	.990	0.991	0.990	.023	.002	.035	−.004
	G4	.969	.984	.993	0.991	0.996	.018	.004	.033	.005
N = 1,000	G1	.982	.984	.999	0.999	1	.011	.002	.017	.003
$ρ = 0$	G2	.979	.989	.993	0.996	0.998	.011	−.004	.016	−.003
	G3	.966	.984	.991	0.995	0.991	.015	.003	.013	.002
	G4	.976	.984	.995	0.995	0.999	.014	.005	.018	.005
N = 1,000	G1	.974	.977	.996	1	0.999	.012	.004	.018	.003
$ρ = 0.5$	G2	.969	.983	.993	0.988	0.998	.014	−.001	.016	−.003
	G3	.964	.988	.986	0.992	0.993	.016	.002	.014	.001
	G4	.973	.988	.989	0.992	0.999	.016	−.003	.021	.002

Note. “RMSE $_{s}$ ” and “RMSE $_{g}$ ” refer to average RMSEs of slipping and guessing parameters. “Bias $_{s}$ ” and “Bias $_{g}$ ” refer to average biases of slipping and guessing parameters. RMSE = root mean squared error; AAR = attribute accuracy rate; PAR = pattern accuracy rate.

Similar with the findings for the recovery of G, the increase of both sample size and test length leads to better recovery of CDM model parameters. Table 6 summarizes the attribute profile accuracy and item parameter estimation from the DINA when $N = 1,000 and J = 40$ . The result is consistent with the simulation studies by Sen and Cohen (2021), where the classification accuracy for DINA model may not be monotonic with increased sample sizes. The recovery of attribute profiles may be jointly affected by sample size, test length, and mastery rate. Table 7 reports the RMSEs of estimated $π$ under the $N = 500 and J = 20$ , $N = 1,000 and J = 20$ and $N = 1,000 and J = 40$ cases. The RMSEs of $\hat{π}$ are smaller for $N = 1,000 and J = 40$ case compared to $N = 500 and J = 20$ cases, which is consistent with the PAR estimations. The average biases of estimated $π$ are less than 0.002 across different sample sizes and test lengths.

Table 6.

Estimated Attribute Profile Accuracy and Item Parameter RMSEs and Biases of the Proposed Method for K = 4 and J = 40 Cases

	Structure	PAR	AAR1	AAR2	AAR3	AAR4	RMSE_s	Bias_s	RMSE_g	Bias_g
$N = 1,000$	G1	.994	.996	.999	1	1	.011	−.003	.016	.001
$ρ = 0$	G2	.993	.995	.998	0.999	0.999	.012	.001	.016	.004
	G3	.991	.995	.998	0.998	0.999	.013	.003	.013	.001
	G4	.993	.996	.999	0.998	0.999	.011	.003	.014	−.001
$N = 1,000$	G1	.985	.988	.996	1	1	.013	−.001	.015	.002
$ρ = 0.5$	G2	.983	.989	.997	0.997	0.999	.013	.002	.013	.004
	G3	.981	.989	.997	0.997	0.998	.014	.002	.014	−.004
	G4	.984	.990	.997	0.997	0.999	.012	.000	.015	.001

Table 7.

RMSEs of Estimated π for K = 4 Cases

	$N = 500 J = 20$		$N = 1,000 J = 20$		$N = 1,000 J = 40$
Structure	$ρ = 0$	$ρ = 0.5$	$ρ = 0$	$ρ = 0.5$	$ρ = 0$	$ρ = 0.5$
G1	.009	.011	.006	.007	.002	.004
G2	.007	.009	.007	.009	.003	.004
G3	.005	.008	.005	.008	.003	.005
G4	.005	.008	.007	.009	.004	.006

Note. The RMSEs are computed by averages of RMSEs of estimated attribute profile proportions across 100 replications. RMSE = root mean squared error.

Recovery of G for K = 7

To simplify the presentation, instead of presenting the $7 \times 7$ adjacency matrix that documents the inclusion probability for each pair of attributes, we presented Figure 7 with the inclusion probabilities annotated by the side of estimated edges when test length equals 20. First, the estimated inclusion probabilities for true edges are similar for $ρ = 0$ and $ρ = 0.5$ , except for the linear structure G ₁. Second, different from what we observe for $K = 4$ , the algorithm estimated a few edges that violate the true hierarchy structure for G ₁, G ₂, and G ₃. For $N = 1,000$ case, the inclusion probabilities are slightly worse than $N = 500$ case. This is because for DINA model, attribute classification accuracy depends on sample size, test length, and attributes together, and the hierarchy structure estimation is affected accordingly. However, the proposed algorithm estimated the existing true edges (and their implied edges) with inclusion probabilities larger than the estimated nonexist edges. In fact, as K gets larger, the sample space for DAGs becomes larger; thus, it is more difficult to detect the true structure for $K = 7$ than for $K = 4$ . Finally, we observe that the linear structures (G ₁ and G ₃) are more difficult to recover than the divergent structure $G_{2}$ and convergent structure G ₄.

Figure 7.

Estimated directed acyclic graph edge probabilities from simulated data for K = 7 and J = 20. Note. The thick edges are the edges from the true structures, and the dashed edges are edges not appearing in the true structures. Implied edges are omitted from the graph.

Table 8 summarizes the overall structure estimation over 100 replications for test length equals 20. Because the inclusion probabilities for $K = 7$ are relatively smaller than $K = 4$ estimations, we set the threshold of inclusion probability for final estimation to be 0.25. As presented in Figure 7, the inclusion probability of true edges is larger than the inclusion probability of nonexisting edges in the simulation results. Therefore, if there are two edges of opposite directions among a pair of attributes, we include the edge with the larger inclusion probability in the final estimation. The proposed method identifies the true structure or a subset of it for $K = 7$ for most replications. The algorithm performs better on the divergent G ₂ and convergent G ₄, compared to linear structures G ₁ and G ₃.

Table 8.

Frequency of Corrected Estimated Structure Among 100 Replications for K = 7 and J = 20 Cases

		$N = 500$			$N = 1,000$
$ρ$	Structure	True	Subset	False	True	Subset	False
$ρ = 0$	G1	82	4	14	72	25	3
	G2	96	4	0	96	4	0
	G3	79	12	9	71	29	0
	G4	91	8	1	75	25	0
$ρ = 0.5$	G1	65	18	17	65	30	5
	G2	90	8	2	88	12	0
	G3	56	18	26	54	33	13
	G4	72	23	5	68	24	8

Note. The threshold of inclusion probability is set to be $0.25$ . “True” refers to estimated structures (including implied edges) equivalent to the true structure, “Subset” refers to estimated structures missing some edges, and “False” refers to estimated structures with extra edges or wrong directions.

Similar to $K = 4$ cases, the performance of attribute hierarchy recovery improves with the increase of both sample size and test length. Figure 8 and Table 9 report the inclusion probability and corrected estimated structures, respectively, for $N = 1,000 and J = 40$ . Compared to $J = 20$ case, there are no incorrect edges; all estimated structures are true or at least a subset of the true structure.

Figure 8.

Table 9.

Frequency of Corrected Estimated Structure Among 100 Replications for K = 7, N = 1,000, and J = 40 Cases

		$N = 1,000$
$ρ$	Structure	True	Subset	False
$ρ = 0$	G1	93	7	0
	G2	94	6	0
	G3	89	11	0
	G4	91	9	0
$ρ = 0.5$	G1	90	10	0
	G2	90	10	0
	G3	85	15	0
	G4	83	17	0

Recovery of CDM parameters for K = 7

Tables 10 and 11 report the overall attribute classification accuracy rate, the estimation accuracy rate for each attribute, and the RMSEs and biases of slipping and guessing parameters. The overall findings are similar to $K = 4$ . However, we note that the overall PAR accuracy rate is smallest on divergent G ₂, which is different from its hierarchy estimation performance. This is due to the fact that G ₂ results in the largest number of permissible attribute patterns among the four structures, which naturally produces lower PAR than others. The AAR of each attribute of G ₂ is comparable with those from all other structures. For $N = 1,000$ case, the attribute classification accuracy is slightly worse than $N = 500$ case, and the RMSEs are smaller for a larger sample size. The average biases of the slipping and guessing parameters are close to 0.

Table 10.

Estimated Attribute Profile Accuracy and Item Parameter RMSEs of the Proposed Method for K = 7 and J = 20 Cases

	Structure	PAR	AAR1	AAR2	AAR3	AAR4	AAR5	AAR6	AAR7	RMSE $_{s}$	Bias $_{s}$	RMSE $_{g}$	Bias $_{g}$
$N = 500$	G1	.927	.983	.996	.996	.990	.995	.969	.988	.018	−.003	.030	.006
$ρ = 0$	G2	.832	.984	.989	.987	.960	.978	.940	.955	.022	−.004	.027	−.008
	G3	.892	.984	.990	.993	.974	.991	.967	.973	.018	.005	.029	.006
	G4	.896	.975	.993	.990	.991	.992	.948	.989	.019	.006	.027	−.008
$N = 500$	G1	.825	.990	.996	.878	.989	.988	.965	.993	.019	.008	.042	−.003
$ρ = 0.5$	G2	.794	.996	.987	.987	.943	.972	.921	.940	.024	.006	.039	.006
	G3	.843	.989	.993	.965	.989	.987	.943	.991	.021	.003	.036	−.006
	G4	.878	.996	.993	.990	.954	.976	.957	.986	.022	.008	.037	−.009
$N = 1,000$	G1	.865	.958	.988	.985	.942	.964	.939	.932	.021	.008	.029	−.006
$ρ = 0$	G2	.854	.987	.987	.976	.933	.961	.928	.941	.022	.007	.026	.007
	G3	.763	.953	.988	.969	.938	.966	.931	.936	.017	.006	.027	−.008
	G4	.848	.969	.988	.980	.971	.976	.912	.930	.019	.009	.026	.006
$N = 1,000$	G1	.788	.967	.988	.991	.977	.989	.910	.928	.019	.006	.032	−.008
$ρ = 0.5$	G2	.781	.960	.984	.990	.973	.975	.927	.929	.022	−.008	.033	.009
	G3	.749	.934	.977	.978	.921	.948	.925	.933	.023	.006	.035	.008
	G4	.769	.979	.987	.979	.942	.961	.922	.943	.022	−.005	.033	.008

Note. “RMSE $_{s}$ ” and “RMSE $_{g}$ ” refer to RMSEs of slipping and guessing parameters. RMSE = root mean squared error; AAR = attribute accuracy rate; PAR = pattern accuracy rate.

Table 11.

Estimated Attribute Profile Accuracy and Item Parameter RMSEs of the Proposed Method for K = 7 and J = 40 Cases

	Structure	PAR	AAR1	AAR2	AAR3	AAR4	AAR5	AAR6	AAR7	RMSE $_{s}$	Bias $_{s}$	RMSE $_{g}$	Bias $_{g}$
N = 1,000	G1	.965	.984	.998	.999	.997	.998	.997	.991	.014	.002	.012	−.004
$ρ = 0$	G2	.909	.992	.996	.998	.988	.994	.963	.969	.013	.003	.015	−.003
	G3	.878	.972	.997	.989	.995	.999	.923	.997	.017	−.006	.016	.007
	G4	.905	.991	.991	.993	.982	.987	.956	.999	.014	.004	.015	.003
N = 1,000	G1	.941	.970	.994	.999	.996	.997	.985	.999	.015	−.007	.019	.006
$ρ = 0.5$	G2	.897	.989	.998	.996	.985	.994	.956	.967	.021	−.005	.025	.008
	G3	.838	.948	.999	.991	.994	.999	.923	.974	.023	.004	.025	.006
	G4	.901	.991	.996	.997	.985	.987	.949	.991	.016	−.006	.021	.005

Note. “RMSE $_{s}$ ” and “RMSE $_{g}$ ” refer to RMSEs of slipping and guessing parameters. RMSE = root mean squared error; AAR = attribute accuracy rate; PAR = pattern accuracy rate.

Table 12 reports the average computation time (in seconds) of the proposed algorithm under each condition. Simulation and estimation work was completed using R by a seamless interfacing established with RcppArmadillo and Rcpp packages (Eddelbuettel et al., 2011). On average, for $K = 4$ , the proposed algorithm took less than 6 minutes for a chain with length $20,000$ , and for $K = 7$ , it takes about 45 minutes to finish a chain with $N = 1,000$ , $J = 40$ and chain length $20,000$ . The computation time increases exponentially in K because of $2^{K}$ possible attribute profile patterns and $K (K - 1)$ possible candidate edges.

Table 12.

Average Computation Time of the Proposed Algorithm Under Each Condition in Seconds

N	K	J	$ρ$	G1	G2	G3	G4
500	4	20	0	156.25	162.86	156.39	156.70
500	4	20	0.5	157.38	162.52	156.92	158.38
1,000	4	20	0	214.47	229.23	229.70	229.75
1,000	4	20	0.5	233.83	240.48	237.79	235.53
1,000	4	40	0	287.44	321.94	307.56	309.55
1,000	4	40	0.5	292.33	324.03	319.99	312.64
500	7	20	0	1,119.36	1,137.94	1,124.62	1,135.18
500	7	20	0.5	1,123.78	1,139.26	1,130.20	1,136.79
1,000	7	20	0	1,944.25	2,090.44	1,967.32	2,075.93
1,000	7	20	0.5	2,082.74	2,131.28	2,065.55	2,187.02
1,000	7	40	0	2,579.00	2,618.58	2,614.96	2,616.03
1,000	7	40	0.5	2,646.88	2,131.28	2,633.92	2,658.55

Note. The computation time is averaged across 100 replications under each condition. The computing is done on a single CPU core with 4-GB memory.

5. Data Application

In this section, we applied the proposed MH-Gibbs algorithm in estimating attribute hierarchy on the Examination for the Certificate of Proficiency in English (ECPE) data. This data set contains 2,922 test takers’ responses to 28 ECPE test items, which has been analyzed by several studies regarding attribute hierarchy (e.g., Ma & Xu, 2021; Templin & Bradshaw, 2014; Wang & Lu, 2021). Three attributes were identified, including (1) morphosynatacic rules, (2) cohesive rules, and (3) lexical rules. The response matrix and Q matrix can be loaded in the R package “CDM” (George et al., 2016). A linear hierarchical structure of the three attributes, $3 \to 2 \to 1$ , was specified or estimated by several previous studies (e.g., Templin & Bradshaw, 2014; Wang & Lu, 2021). In this analysis, we applied the developed algorithm to learn the attribute hierarchy structure from the data using both DINA and GDINA model. For each CDM model, a total of five chains with a random starting point were first ran to check the convergence using the maximum value of Gelman–Ruin proportional scale reduction factor (Gelman & Rubin, 1992). Results indicate the algorithm converged after 3,000 iterations. Thus, we summarize the result based on a single chain of the last 5,000 iterations of the 10,000 iterations. The estimated adjacency matrix from both DINA and GDINA was very consistent, and here, we presented the result based on the DINA model in Table 13. We can conclude that Attribute 3 is the prerequisite of Attribute 2 as the (3, 2) element is nonzero (0.584). Thus, our algorithm estimated a subgraph of the expert specified attribute hierarchy, which is consistent with the finding from the simulation study that our algorithm is able to at least recover a subgraph of the true structure. This demonstrates that our estimated structure can serve as the starting structure for cognitive experts if no prior information regarding the attribute hierarchy exists. The proposed method did not identify the relationship of $α_{c}$ compared to the previous studies. We found the performance of the proposed algorithm could potentially be affected by the prespecified Q matrix, where the three skills were loaded unbalanced among the test items. The missing edge allows two more permissible patterns (1,0,0) and (1,0,1), and the corresponding ideal response pattern of these two attribute profiles only differs from other permissible patterns by 5 items under the currently prespecified unstructured Q matrix. Further improvement of the proposed algorithm may require joint modeling of attribute structure and corresponding structured Q matrix.

Table 13.

The Estimated Adjacency Matrix for Examination for the Certificate of Proficiency in English Data

	1	2	3
1	0	0.010	0
2	0	0	0
3	0	0.584	0

6. Discussion

Attribute hierarchy plays an important role in applying CDMs for designing efficient cognitive diagnostic assessments, yet there are limited statistical tools to directly learn the latent hierarchy structure among attributes. In this study, we formulated the estimation of attribute hierarchy underlying CDMs with a Bayesian framework and developed an efficient MCMC algorithm to estimate the hierarchy structure together with CDM model parameters directly. The exploration of sampling space of hierarchy structure uses a random walk on transitive reduction of acyclic directed graphs to effectively search possible hierarchy states. Using the DINA model as an example, our simulation study demonstrates that the proposed method can fully recover or estimate at least a subgraph of the underlying structure across various conditions. In addition, we applied our proposed algorithm in estimating the hierarchy using the ECPE data. Results indicate that our algorithm is able to recover a subgraph of the prespecified attribute hierarchy. This indicates the potential of using our developed algorithm in learning attribute hierarchy from data, and the results can be used to validate the existing attribute hierarchy specified by content experts.

One finding from the simulation study is about the impact of the true attribute hierarchical structure on our proposed estimation algorithm. We note that the linear structure is more difficult to estimate correctly than nonlinear structures. This may be due to the construction of the M–H step in the proposed algorithm, where we proposed to only add or remove one edge at a time. Other stochastic construction of graphs, such as exchanging edges or reverting direction, can be considered in the algorithm to further expedite the exploration of the sample space of DAGs in the future.

Another finding from the simulation study is that the inclusion probabilities calculated from sampled posterior adjacency matrices are often smaller than $0.5$ , which indicates that the acceptance ratio of adding an edge to the free structure is not very high. This is because currently we used a Dirichlet prior for the proportions of attribute profiles based on the nonhierarchical structure, and the hyperparameter in simulation and data application is $δ = (1, \dots {,1)}^{Τ}$ . When the current G is the nonhierarchical structure (all zeros) and the newly proposed $G^{*}$ in Algorithm 2 Line 6 has one edge, $\frac{p (Y | π, G^{*})}{p (Y | π, G)} = 1$ , then $\frac{p (π | G^{*})}{p (π | G)}$ dominates the acceptance ratio. Therefore, in practice, one can use an unbalanced $δ_{0}$ to deliberately make the acceptance ratio in favor of certain dependent structures. This can further improve the inclusion probabilities for potential edges.

The proposed method gives an exploratory discovery of potential hierarchical structures. The estimated adjacency matrix with inclusion probability of potential edges provides more information on the possible prerequisite relationships among attributes than a strict estimated structure does. However, if the researchers look for further interpretation on the attribute mastery, an explicit structure is needed. The choice of a threshold provides guidance on how to decide on explicit structures. The selection of threshold values to include the edge is subjective in this study. A potential strategy is to consider the overall possible DAGs of size K. The exact enumeration given by Robinson (1977) shows the total number of possible DAGs with K nodes is given by the recurring function $G (K)$ :

G (K) = \sum_{k = 1}^{K} {(- 1)}^{k - 1} (\begin{matrix} K \\ k \end{matrix}) 2^{k (K - k)} G (K - k) .

Because the proposed algorithm considers a uniform prior, then the edges (off-diagonal elements in G) have an equal probability of appearing in the sample space of DAGs, which is upper bounded by $\frac{5^{K - 2} G (K - 2)}{G (K)}$ . Because $G (K)$ is superexponential in K, therefore one can consider an approximate value of $e^{- 2}$ as the threshold. Alternatively, a subjective threshold depending on posterior samples such as a certain percentile of edge inclusion probability could also be considered.

Even though the current study mainly used the DINA model in the simulation study, the algorithm can be generally applicable to other CDM specifications as long as the posteriors can be derived explicitly. The Appendix showed another simple simulation of the proposed method on a compensatory main-effect restricted latent class model. In addition, one limitation of the current study is that we assume Q is known and unstructured. In the future, we can estimate Q together with other CDM model parameters and the underlying attribute hierarchy. Our current model framework assumes that the item parameter matrix B fully implies the structure of Q; thus, the estimation of structured Q can be easily incorporated into the existing algorithm through the estimated sparsity of B. Furthermore, the proposed algorithm uses a Dirichlet prior on attribute pattern proportions $π$ and utilizes the aggregation property of Dirichlet distribution. An alternative formulation using Polya urn model for attribute profiles can be further investigated to collapse the parameter $π$ . Finally, future research can potentially extend the proposed algorithm to longitudinal CDMs, where the transition probabilities among different attribute patterns are also constrained by G. Under specific model assumptions such as the first-order Markov model, the proposed algorithm is also feasible with explicit derivation of the transition matrix analogous to the construction of $π$ .

Footnotes

Appendix

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research and/or authorship of this article: This research was supported by National Science Foundation under award SES-2051198.

ORCID iD

Yinghan Chen

References

Aho

A. V.

Garey

M. R.

Ullman

J. D.

(1972). The transitive reduction of a directed graph. SIAM Journal on Computing, 1(2), 131–137.

Brooks

S. P.

Gelman

(1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4), 434–455.

Chen

Culpepper

Liang

(2020). A sparse latent class model for cognitive diagnosis. Psychometrika, 85(1), 121–153.

Chen

Liu

Ying

(2018). Recommendation system for adaptive learning. Applied Psychological Measurement, 42(1), 24–41.

Chen

Liu

Ying

(2015). Statistical analysis of Q-matrix based diagnostic classification models. Journal of the American Statistical Association, 110(510), 850–866. https://doi.org/10.1080/01621459.2014.934827

Chickering

D. M.

(2002). Learning equivalence classes of Bayesian-network structures. Journal of Machine Learning Research, 2(1), 445–498.

Colombo

Maathuis

M. H.

(2014). Order-independent constraint-based causal structure learning. The Journal of Machine Learning Research, 15(1), 3741–3782.

Colombo

Maathuis

M. H.

Kalisch

Richardson

T. S.

(2012). Learning high-dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, 40(1), 294–321.

Culpepper

S. A.

(2015). Bayesian estimation of the DINA model with Gibbs sampling. Journal of Educational and Behavioral Statistics, 40(5), 454–476.

10.

de la Torre

(2011). The generalized DINA model framework. Psychometrika, 76(2), 179–199. https://doi.org/10.1007/s11336-011-9207-7

11.

de la Torre

Douglas

J. A.

(2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69(3), 333–353. https://doi.org/10.1007/bf02295640

12.

de la Torre

Douglas

J. A.

(2008). Model evaluation and multiple strategies in cognitive diagnosis: An analysis of fraction subtraction data. Psychometrika, 73(4), 595–624.

13.

de la Torre

Hong

Deng

(2010). Factors affecting the item parameter estimation and classification accuracy of the DINA model. Journal of Educational Measurement, 47(2), 227–249.

14.

Eddelbuettel

François

Allaire

Chambers

Bates

Ushey

(2011). Rcpp: Seamless R and C++ integration. Journal of Statistical Software, 40(8), 1–18.

15.

Gelman

Rubin

D. B.

(1992). Inference from iterative simulation using multiple sequences. Statistical Science, 457–472. https://doi.org/10.1214/ss/1177011136

16.

George

A. C.

Robitzsch

Kiefer

Groß

Ünlü

(2016). The R package CDM for cognitive diagnosis models. Journal of Statistical Software, 74, 1–24.

17.

Leighton

J. P.

Gierl

M. J.

(2007). Why cognitive diagnostic assessment. In Leighton

J. P.

Gierl

M. J.

(Eds.), Cognitive diagnostic assessment for education: Theory and applications (pp. 3–18). Cambridge University Press. https://doi.org/10.1017/cbo9780511611186.001

18.

Leighton

J.P.

Gierl

Hunka

(2004). The attribute hierarchy model: An approach for integrating cognitive theory with assessment practice. Journal of Educational Measurement, 41(3), 205–236.

19.

Ouyang

(2022). Learning latent and hierarchical structures in cognitive diagnosis models. Psychometrika, 88, 175–207.

20.

(2021). Hypothesis testing for hierarchical structures in cognitive diagnosis models. Journal of Data Science, 20(3), 279–302.

21.

Robinson

R. W.

(1977). Counting unlabeled acyclic digraphs. In Combinatorial mathematics v (pp. 28–43). Springer. https://doi.org/10.1007/BFb0069178

22.

Sen

Cohen

A. S.

(2021). Sample size requirements for applying diagnostic classification models. Frontiers in Psychology, 11, 621251.

23.

Sorrel

M. A.

Olea

Abad

F. J.

de la Torre

Aguado

Lievens

(2016). Validity and reliability of situational judgement test scores: A new approach based on cognitive diagnosis models. Organizational Research Methods, 19(3), 506–532.

24.

Tatsuoka

(2002). Data analytic methods for latent partially ordered classification models. Journal of the Royal Statistical Society: Series C (Applied Statistics), 51(3), 337–350.

25.

Tatsuoka

K. K.

(1984). Analysis of errors in fraction addition and subtraction problems. Computer-Based Education Research Laboratory.

26.

Tatsuoka

K. K.

(1985). A probabilistic model for diagnosing misconceptions by the pattern classification approach. Journal of Educational and Behavioral Statistics, 10(1), 55–73.

27.

Tatsuoka

K. K.

(2009). Cognitive assessment: An introduction to the rule space method. Routledge.

28.

Templin

J. L.

Bradshaw

(2014). Hierarchical diagnostic classification models: A family of models for estimating and testing attribute hierarchies. Psychometrika, 79(2), 317–339.

29.

Wang

Cai

Douglas

Chang

H.-H.

(2019). Cognitive diagnostic models with attribute hierarchies: Model estimation with a restricted Q-matrix design. Applied Psychological Measurement, 43(4), 255–271.

30.

von Davier

(2008). A general diagnostic model applied to language testing data. British Journal of Mathematical and Statistical Psychology, 61(2), 287–307.

31.

Wang

Gierl

M. J.

(2011). Using the attribute hierarchy method to make diagnostic inferences about examinees’ cognitive skills in critical reading. Journal of Educational Measurement, 48(2), 165–187.

32.

Wang

(2021). Learning attribute hierarchies from data: Two exploratory approaches. Journal of Educational and Behavioral Statistics, 46(1), 58–84.