Bayesian Analysis Methods for Two-Level Diagnosis Classification Models

Abstract

Understanding whether or not different types of students master various attributes can aid future learning remediation. In this study, two-level diagnostic classification models (DCMs) were developed to represent the probabilistic relationship between external latent classes and attribute mastery patterns. Furthermore, variational Bayesian (VB) inference and Gibbs sampling Markov chain Monte Carlo methods were developed for parameter estimation of the two-level DCMs. The results of a parameter recovery simulation study show that both techniques appropriately recovered the true parameters; Gibbs sampling in particular was slightly more accurate than VB, whereas VB performed estimation much faster than Gibbs sampling. The two-level DCMs with the proposed Bayesian estimation methods were further applied to fourth-grade data obtained from the Trends in International Mathematics and Science Study 2007 and indicated that mathematical activities in the classroom could be organized into four latent classes, with each latent class connected to different attribute mastery patterns. This information can be employed in educational intervention to focus on specific latent classes and elucidate attribute patterns.

Keywords

diagnostic classification models latent class analysis variational Bayesian inference Gibbs sampling algorithm

1. Introduction

Diagnosis classification models (DCMs; e.g., Rupp et al., 2010; von Davier & Lee, 2019) are educational measurement models that enable learners and teachers to focus on cognitively weak points and improve the individual learning status of the students. Well-known general models such as the log-linear cognitive diagnostic model (Henson et al., 2009), generalized deterministic input noisy and-gate model (de la Torre, 2011), and general diagnostic model (von Davier, 2008) have been applied to real data such as those obtained by the Programme for International Student Assessment (Chen & de la Torre, 2014) or Trends in International Mathematics and Science Study (TIMSS; e.g., Yamaguchi & Okada, 2018). These analyses have revealed what types of fine cognitive elements, called attributes in the DCM context, are crucial for answering test items. However, DCMs typically use only response data for test items, and their item response and structural parameters represent the assumptions of diagnostic tests. Furthermore, general DCMs do not reveal how attributes are connected with external variables.

One type of extension of DCMs that includes external variables is called explanatory DCM (e.g., Park & Lee, 2014; Park et al., 2018). Explanatory DCMs assume that the attribute mastery probabilities or item response probabilities are predicted with external variables. Through them, one can understand how prediction can be performed for a mathematical attribute, such as the addition or subtraction of fractions, by using the learning time of an individual or their attitude or motivation to learn the attribute. Explanatory DCMs can explicitly model attribute mastery and correct item responses affected by different factors from latent attributes. Park et al. (2018) investigated the effects of mathematics confidence (as measured by items such as “I usually do well in mathematics”), affect in mathematics (as measured by items such as “I enjoy learning mathematics”), and calculator ownership status (1 = having, 0 = not having) as predictors. The results summarized in Table 3 by Park et al. (2018) indicated that owning a calculator and feeling more confident decreased the probability of mastering attributes such as “whole numbers” or “dimensions and locations.” However, an increase in the score of affect in mathematics indicated an increase in the probabilities of mastering attributes. In another related model, Zhan et al. (2018) proposed DCMs incorporating response time, which is another source of information about the latent abilities of individuals. Zhan et al. (2018) developed a hierarchical DCM with higher order latent traits that are associated with response time. Other types of external information were employed by Zhan et al. (2022), who proposed a DCM incorporating not only response times but also biometrics data such as visual fixation counts. These models provide flexible frameworks to deal with additional sources of information to diagnose individual attribute mastery status.

However, explanatory DCMs do not consider the connections between the attribute mastery pattern level and explanatory variables; they instead focus on each attribute level model. This implies that the explanatory DCMs proposed by Park et al. (2018) assume conditional independence among attributes given predictors. Therefore, the attribute mastery pattern probabilities are represented as the product of conditional probability of attributes given covariates. The assumption is that the attribute masteries are independent of each other when the predictors are fixed; however, this assumption is too strong because it is not always possible to include covariates to completely explain the relationships among attributes. Instead, modeling attribute mastery pattern level can moderate the assumption.

Furthermore, the explanatory DCMs mainly use continuous variables as covariates rather than categorical latent variables known as latent classes (e.g., Collins & Lanza, 2009; Hagenaars & McCutcheon, 2002; Lazarsfeld & Henry, 1969; White & Murphy, 2014). Both continuous observed and latent variables are generally employed as covariates (e.g., De Boeck, 2004). However, latent classes that connect with attribute mastery patterns can be more useful in group-level educational remedies in a classroom setting.

For example, in math classes, “activeness” forms latent classes, which are active and nonactive commitments determined by the various actual activities in a math class such as involving group discussion or speaking up in class. Here, the nonactive commitment class may be strongly connected to low attribute mastery patterns, while the active commitment class is related to high mastery patterns. In this case, teachers can consider methods to encourage students to be involved in class activities to improve their mathematical attribute mastery. In essence, one latent class layer explains the attribute mastery patterns and the latent class can provide hints to educational remedy. In other words, the attribute mastery patterns and latent classes consist of a hierarchical relationship, and these latent classes can reveal the types of individuals that tend to belong to specific attribute mastery patterns.

DCMs are obtained by constraining general latent class models (Rupp & Templin, 2008), and attribute mastery patterns are latent classes. Therefore, assuming an additional latent class and assessing the relationships between the latent classes and attribute mastery patterns can be considered a special case of a two-level latent class model (Miyazaki et al., 2007). A two-level latent class can directly model the strength of connection between two latent classes as probability parameters. This idea can be applied to the DCM context, and external variables can be considered differently than explanatory DCMs.

In this research, two-level DCMs whose attribute mastery patterns were connected to exogenous latent classes were developed. Latent class–related models sometimes exhibit poor parameter estimation because of the sparseness of cells. To this end, utilizing a Bayesian prior can stabilize the estimation in latent class–related models (Collins & Lanza, 2009, p. 171). Therefore, we developed two Bayesian estimation methods for the parameter estimation of two-level DCMs: the variational Bayesian (VB) inference and Gibbs sampling methods. The VB inference method (e.g., Ormerod & Wand, 2010) is faster than the Markov chain Monte Carlo (MCMC) approach (e.g., Brooks et al., 2011) but provides parameter estimates similar to a fully Bayesian estimation method. The VB method has been applied to test models such as general DCMs (e.g., Yamaguchi & Okada, 2020a). Generally, the MCMC method is computationally heavy; however, the Gibbs sampling method is a more effective MCMC method than the Metropolis–Hastings type MCMC and can approximate posteriors more precisely than the VB method. Thus, to examine the quality of parameter estimates, we compared these Bayesian methods through a simulation study.

The contribution of this study includes not only developing a DCM that enables us to represent a connection between an attribute mastery pattern and external latent classes but also developing stable Bayesian estimation methods. An expectation–maximization (EM) algorithm to achieve classical maximum-a-posteriori (MAP) estimation, which is a common point-estimation-based Bayesian method, was also derived and is presented in Online Appendix A. The remainder of this article is structured as follows: Section 2 introduces the formulation of two-level DCMs and derivations of the aforementioned Bayesian estimation methods. Section 3 presents a simulation study to check the parameter recovery of the proposed estimation methods. Section 4 provides a real data demonstration of two-level DCMs with TIMSS data. Finally, Section 5 summarizes the conclusions, provides a discussion, and mentions future research possibilities.

2. Formulation of Two-Level DCMs and Bayesian Estimation Procedure

2.1. Elements of Two-Level DCMs

The core idea of two-level DCMs is to combine DCMs and latent class models using relationship probabilities that connect both measurement models. Figure 1 provides a conceptual representation of two-level DCMs, which have three elements: diagnosis measurement, latent class measurement, and structural models. The essential element in two-level DCMs is the structural model among the exogenous latent class and attribute mastery patterns. This structural model represents the strength of the connection between the latent classes and attribute mastery patterns. These models are introduced hereafter. In addition, to derive Bayesian estimation methods, priors for the model parameters are specified, and joint distributions of the observed variables, latent variables, and parameters are defined.

Figure 1.

Conceptual representation of the two-level diagnostic classification model.

For the diagnostic measurement, we employed the latent class formulation and notations utilized in Yamaguchi and Okada (2020a) and Yamaguchi and Templin (2022a, 2022b). Attribute mastery patterns $α = {[α_{1}, \dots, α_{k}, \dots, α_{K}]}^{⊤}$ are defined as permutations of the mastery status of each attribute. If the $k (= 1, \dots, K)$ th attribute $α_{k}$ is mastered, then $α_{k}$ is 1; otherwise, it is 0. In this study, all possible mastery patterns were considered; thus, $2^{K}$ possible permutations could be assumed. Therefore, let $l (= 1, \dots, 2^{K})$ be a subscript of latent attribute mastery patterns, and $α_{l}$ be the lth attribute mastery pattern, which takes element $α_{l k} \in {0, 1}$ , that is, the kth attribute mastery status of the lth attribute. $α_{l k} = 1 (= 0)$ indicates the mastery (nonmastery) of the attribute. In addition, one-of- $2^{K}$ coding of the latent attribute mastery patterns of individual $i (= 1, \dots, I)$ is $z_{i} = {[z_{i 1}, \dots, z_{i l}, \dots, z_{i 2^{K}}]}^{⊤}$ , whose element $z_{i l} \in {0, 1}$ indicates that individual i belongs to the lth attribute mastery pattern if $z_{i l} = 1$ ; otherwise, it is 0. Therefore, _i satisfies $\sum_{l} z_{i l} = 1$ . By using attribute mastery indicator z _i , DCMs can be clearly expressed as latent class (mixture) models.

Next, we define the probability of the $j (= 1, \dots, J)$ th item response probability of individual i given attribute mastery pattern indicator z _i . To define the correct item response probabilities of DCMs, it is necessary to introduce a Q-matrix and constraint matrices. The Q-matrix represents which attributes are measured by which item; thus, an element is $q_{j k} \in {0, 1}$ , which is 1 if item j measures the kth attribute and 0 otherwise. In general DCMs, this Q-matrix generates item-specific patterns using the attributes measured by the item. For example, if the jth item measures the first and second items and its q-vector is $q_{j} = [1, 1, 0]$ , then the item can separate four possible attribute mastery patterns, which are named according to the specific pattern $α_{j h}^{*}$ related to the first and second attributes: $α_{j 1}^{*} = [0, 0, *]$ , $α_{j 2}^{*} = [0, 1, *], α_{j 3}^{*} = [1, 0, *]$ , and $α_{j 4}^{*} = [1, 1, *]$ , where $h (= 1, \dots, 2^{\sum_{k} q_{j k}})$ represents the item-specific pattern. The asterisks mean that the third attribute is not considered in these items. The constraint matrix G_j has elements $g_{j h l} \in {0, 1}$ , where $g_{j h l}$ is 1 if all the nonasterisk elements of $α_{j h}^{*}$ are equal to the corresponding elements of $α_{l}$ and 0 otherwise. This G_j matrix indicates that several attribute mastery patterns are reduced to one of the item-specific patterns. Detailed definitions and examples are shown in Yamaguchi and Okada (2020a).

Using G_j and z _i , the correct response probability of an item-specific pattern is represented as $θ_{j h} = P (X_{i j} = 1 | \sum_{l} z_{i l} g_{j h l} = 1)$ , and the conditional likelihood of individual i for $θ_{j} = {θ_{j h}}_{h = 1}^{2^{\sum_{k} q_{j k}}}$ given latent attribute mastery pattern indicator z_i is

P (X_{i j} = x_{i j} | θ_{j}, z_{i}, G_{j}) = \prod_{h = 1}^{2^{\sum_{k} q_{j k}}} \prod_{l = 1}^{2^{K}} {[{θ_{j h}^{x_{i j}} {(1 - θ_{j h})}^{1 - x_{i j}}}^{g_{j h l}}]}^{z_{i l}} .

Assuming conditional independence and exchangeability, the conditional probability of $X = {x_{i j}}_{i = 1, j = 1}^{I, J}$ given $Θ = {θ_{j}}_{j = 1}^{J}$ , $Z = {z_{i}}_{i = 1}^{I}$ , and $G = {G_{j}}_{j = 1}^{J}$ is

P (X|Θ, Z, G) = \prod_{i = 1}^{I} \prod_{j = 1}^{J} P (X_{i j} = x_{i j} | θ_{j}, z_{i}, G_{j}) .

This is the measurement model of the diagnostic classification approach.

The measurement model of the latent class model assumes unobserved latent groups underlying several categorical indicators (e.g., White & Murphy, 2014). The random categorical indicator $Y_{i p}$ takes a binary value, 0 or 1, for simplicity, and its realization is $y_{i p}$ , where $p (= 1, \dots, P)$ is the number of variables. In addition, the latent class indicator vector $ω_{i} = {[ω_{i 1}, \dots, ω_{i c}, \dots, ω_{i C}]}^{⊤}$ is introduced to represent that individual i belongs to latent class $c (= 1, \dots, C)$ if $ω_{i c} = 1$ and 0 if not. Under the usual latent class model assumptions, such as local independence, and by representing the item response probability parameter as $λ_{p c} = P (Y_{i p} = 1 {|ω}_{i c} = 1)$ , the conditional probability is

P (Y_{i p} = y_{i p} {|λ}_{p c}, ω_{i}) = \prod_{c = 1}^{C} {λ_{p c}^{y_{i p}} {(1 - λ_{p c})}^{1 - y_{i p}}}^{ω_{i c}} .

The conditional probability of $Y = {y_{i p}}_{i = 1, p = 1}^{I, P}$ given $Λ = {λ_{p c}}_{p = 1, c = 1}^{P, C}$ and $Ω = {ω_{i}}_{i = 1}^{I}$ is

P (Y|Λ, Ω) = \prod_{i = 1}^{I} \prod_{p = 1}^{P} P (Y_{i p} = y_{i p} {|λ}_{p c}, ω_{i}) .

The final part is the structural model between the individual attribute mastery pattern indicator z _i and latent class indicator $ω_{i}$ can be expressed as follows:

P (z_{i} | ω_{i}, τ_{c}) = \prod_{l = 1}^{2^{K}} τ_{c l}^{ω_{i c} z_{i l}},

where $τ_{c l} = P (z_{i l} = 1 | ω_{i c} = 1)$ , implying the probability of individual i belonging to the lth attribute mastery pattern when they are the member of the cth latent class, and $τ_{c}$ is a parameter set of relationship probability ${τ_{c l}}_{l = 1}^{2^{K}}$ . $τ_{c l}$ is the strength of the relationship between the latent class and attribute mastery pattern. $τ_{c l}$ satisfies $\sum_{l} τ_{c l} = 1$ and can be interpreted as the strength of the connection between latent class c and attribute mastery pattern l. $τ_{c l}$ works similarly to the transition probability in latent transition analysis. Exchangeability of individuals leads to the following joint distribution:

P (Z|Ω, Τ) = \prod_{i = 1}^{I} \prod_{c = 1}^{C} P (z_{i} | ω_{i}, τ_{c}),

where $Τ = {τ_{c}}_{c = 1}^{C}$ , and we refer to $Τ$ as a relationship matrix. The distribution of latent class indicator $ω_{i}$ should be specified, and it is a categorical or one-trial multinomial distribution with parameter $π = {π_{c}}_{c = 1}^{C}$ , which satisfies $\sum_{c} π_{c} = 1$ :

P (Ω| π) = \prod_{i = 1}^{I} \prod_{c = 1}^{C} π_{c}^{ω_{i c}} .

Furthermore, priors should be set for all parameters. Considering conditional conjugacy, the correct item response probability parameter $θ_{j h}$ and conditional-response probability $λ_{p c}$ have beta distributions with positive-valued parameters $a_{j h}^{0}$ , $b_{j h}^{0}$ , $c_{p c}^{0}$ , and $d_{p c}^{0}$ :

P (θ_{j h} | a_{j h}^{0}, b_{j h}^{0}) = \frac{1}{B e t a (a_{j h}^{0}, b_{j h}^{0})} θ_{j h}^{a_{j h}^{0} - 1} {(1 - θ_{j h})}^{b_{j h}^{0} - 1}, \forall j, h,

P (λ_{p c} | c_{p c}^{0}, d_{p c}^{0}) = \frac{1}{B e t a (c_{p c}^{0}, d_{p c}^{0})} λ_{p c}^{c_{p c}^{0} - 1} {(1 - λ_{p c})}^{d_{p c}^{0} - 1}, \forall p, c,

where $B e t a (x, y)$ is a beta function defined as $Γ (x) Γ (y) / Γ (x + y)$ , and $Γ (\cdot)$ is a gamma function. Assuming independence of the prior distributions and introducing hyperparameter sets $A^{0} = {a_{j h}^{0}}_{j = 1, h = 1}^{J, 2^{\sum_{k} q_{j k}}}$ , $B^{0} = {b_{j h}^{0}}_{j = 1, h = 1}^{J, 2^{\sum_{k} q_{j k}}}$ , $C^{0} = {c_{j h}^{0}}_{p = 1, c = 1}^{P, C}$ , and $D^{0} = {d_{j h}^{0}}_{p = 1, c = 1}^{P, C}$ , the joint distributions of priors are $P ({Θ|A}^{0}, B^{0}) = \prod_{j} \prod_{h} P (θ_{j h} | a_{j h}^{0}, b_{j h}^{0})$ and $P ({Λ|C}^{0}, D^{0}) = \prod_{p} \prod_{c} P (λ_{p c} | c_{p c}^{0}, d_{p c}^{0})$ . In addition, $τ_{c}$ and $π$ are the vectors with probability elements; thus, it is natural to assume Dirichlet distributions with positive-valued element parameter vectors $γ_{c}^{0} = {[γ_{c 1}^{0}, \dots, γ_{c l}^{0}, \dots, γ_{c 2^{K}}^{0}]}^{⊤}$ and $δ^{0} = {[δ_{1}^{0}, \dots, δ_{c}^{0}, \dots, δ_{C}^{0}]}^{⊤}$ ,

P (τ_{c} | γ_{c}^{0}) = \frac{1}{B e t a (γ_{c}^{0})} \prod_{l = 1}^{2^{K}} τ_{c l}^{γ_{c l}^{0} - 1}, \forall c,

P (π| δ^{0}) = \frac{1}{B e t a (δ^{0})} \prod_{c = 1}^{C} π_{c}^{δ_{c}^{0} - 1},

where $B e t a (x = [x_{1}, \dots, x_{m}, \dots, x_{M}]), x_{m} > 0, \forall m$ is a multivariate version of the beta function defined as $\prod_{m = 1}^{M} Γ (x_{m}) / Γ (\sum_{m = 1}^{M} x_{m})$ , which includes the usual beta function defined above. Thus, we use the same symbols for both beta functions and distinguish them by their arguments. Then, the joint distribution of $Τ = {τ_{c}}_{c = 1}^{C}$ is $P ({Τ|Γ}^{0} = {γ_{c}^{0}}_{c = 1}^{C}) = \prod_{c} P (τ_{c} | γ_{c}^{0})$ .

These settings provide the following joint distribution of observed and latent variables and parameters:

P (X, Y, Z, Ω, Θ, Λ, Τ, π {|A}^{0}, B^{0}, C^{0}, D^{0}, Γ^{0}, δ^{0}, G)

= P (X|Θ, Z, G) P ({Θ|A}^{0}, B^{0}) P (Y|Λ, Ω) P ({Λ|C}^{0}, D^{0}) P (Z|Ω, Τ) P ({Τ|Γ}^{0}) P (Ω| π) P (π | δ^{0}) .

This joint distribution was employed to derive the VB inference and Gibbs sampling algorithms, as described in Sections 2.2 and 2.3, respectively.

2.2. VB Inference Algorithm

2.2.1. Basic principle of VB inference

The VB method (Bishop, 2006; Blei et al., 2017; Jeon et al., 2017; Nakajima et al., 2019; Yamaguchi & Okada, 2020a, 2020b) is a popular estimation method in machine learning and has also been employed in psychometrics. Maximum likelihood estimation is a gold standard parameter estimation method in psychometrics but is not always appropriate in situations with small sample sizes or complex models. Bayesian estimation with MCMC is also employed to obtain approximate posteriors but is computationally intractable in situations with large sample sizes or large numbers of parameters. The VB method is computationally tractable but provides Bayesian estimates and does so in a short estimation time. We briefly explain the concept of the VB method and provide a well-established formula to calculate the variational posteriors based on Bishop (2006).

In VB estimation, the variational posterior $q (\cdot)$ for a set of parameters and latent variables $ϑ$ is decomposed into independent variational distribution $q_{m} (\cdot)$ taking mutually exclusive parameter subsets $ϑ_{1}, \dots, ϑ_{m}, \dots, ϑ_{M} (m = 1, \dots, M)$ as arguments:

q (ϑ) = \prod_{m = 1}^{M} q_{m} (ϑ_{m}) .

This decomposition is called mean-field approximation. In addition, let $P (ϑ, X)$ be a joint distribution of parameters, latent variables, and data $X$ . Then, the general formula of variational posterior $q_{m} (ϑ_{m})$ becomes

q_{m} (ϑ_{m}) \propto exp (E_{q_{m^{'} \neq m}} (log P (ϑ, X))),

where $E_{q_{m^{'} \neq m}} (x)$ is $\int \dots \int x \prod_{m^{'} \neq m} q_{m^{'}} (ϑ_{m^{'}}) d ϑ_{1} \dots d ϑ_{m - 1} d ϑ_{m + 1} \dots d ϑ_{M}$ . The VB posteriors depend on each other; therefore, an iterative update procedure is necessary to obtain optimal VB posteriors.

2.2.2. Construction of VB inference algorithm on two-level DCMs

In two-level DCMs, the variational posterior can be expressed as

q (Θ, Λ, Τ, π, Z, Ω) = q (Z, Ω) q (Θ, Λ, Τ, π),

= (\prod_{i = 1}^{I} q (z_{i}, ω_{i})) (\prod_{j = 1}^{J} \prod_{h = 1}^{2^{\sum_{k} q_{j k}}} q (θ_{h j})) (\prod_{p = 1}^{P} \prod_{c = 1}^{C} q (λ_{p c})) (\prod_{c = 1}^{C} q (τ_{c})) q (π) .

The second line in Equation 15 can be naturally derived from the independence between latent variables ${Z, Ω}$ and model parameters ${Θ, Λ, Τ, π}$ . This means that the minimal independence assumption between the latent variables and other parameters leads to further decomposition due to the model structure. The decomposition will provide a simple well-known distribution form for the variational posteriors. The actual distribution forms and VB algorithm, which is also called the variational EM algorithm, will be explained in the next section.

In addition to the variational posteriors, the evidence lower bound (ELBO; i.e., lower bound of log-likelihood) is required to evaluate model appropriateness and a stopping rule for the estimation algorithm. The ELBO can be calculated to assess convergence and is defined as

ELBO = \sum_{Z} \sum_{Ω} \int \int \int \int q (Θ, Λ, Τ, π, Z, Ω) log \frac{P (X, Y, Z, Ω, Θ, Λ, Τ, π)}{q (Θ, Λ, Τ, π, Z, Ω)} d Θ d Λ d Τ d π .

Here, the joint distribution $P (X, Y, Z, Ω, Θ, Λ, Τ, π)$ is Equation 12 with the conditioning variables omitted for simplicity and the variational distribution $q (Θ, Λ, Τ, π, Z, Ω)$ decomposed as in Equation 15. The detailed expression of ELBO is provided in Online Appendix B.

2.2.3. Variational posteriors and VB algorithm

To derive the variational posterior, we borrowed the results of Yamaguchi and Okada (2020a); then, the variational posterior of the correct item response probability parameter becomes a beta distribution again:

q (θ_{j h}) = \frac{1}{B e t a (a_{j h}^{*}, b_{j h}^{*})} θ_{j h}^{a_{j h}^{*}} {(1 - θ_{j h})}^{b_{j h}^{*}}, \forall j, h,

where the parameters are

{\begin{matrix} a_{j h}^{*} = \sum_{i = 1}^{I} \sum_{l = 1}^{2^{K}} E (z_{i l}) g_{j h l} x_{i j} + a_{j h}^{0}, \\ b_{j h}^{*} = \sum_{i = 1}^{I} \sum_{l = 1}^{2^{K}} E (z_{i l}) g_{j h l} (1 - x_{i j}) + b_{j h}^{0} . \end{matrix}

Hereinafter, we take expectation with respect to the VB posteriors. The variational posterior of $λ_{p c}$ is also a beta distribution

q (λ_{p c}) = \frac{1}{B e t a (c_{p c}^{*}, d_{p c}^{*})} λ_{p c}^{c_{p c}^{*} - 1} {(1 - λ_{p c})}^{d_{p c}^{*} - 1}, \forall j, h,

where the parameters are

{\begin{matrix} c_{p c}^{*} = \sum_{i = 1}^{I} E (ω_{i c}) y_{i p} + c_{p c}^{0}, \\ d_{p c}^{*} = \sum_{i = 1}^{I} E (ω_{i c}) (1 - y_{i p}) + d_{p c}^{0} . \end{matrix}

The variational posteriors of $τ_{c}$ and $π$ are both Dirichlet distributions, and their formulas are

q (τ_{c}) = \frac{1}{B e t a (γ_{c}^{*})} \prod_{l = 1}^{2^{K}} τ_{c l}^{γ_{c l}^{*} - 1}, \forall c,

q (π) = \frac{1}{B e t a (δ^{*})} \prod_{c = 1}^{C} π_{c}^{δ_{c}^{*} - 1},

where $γ_{c l}^{*}$ and $δ_{c}^{*}$ are

γ_{c l}^{*} = \sum_{i = 1}^{I} E (ω_{i c} z_{i l}) + γ_{c l}^{0},

δ_{c}^{*} = \sum_{i = 1}^{I} E (ω_{i c}) + δ_{c}^{0} .

Variational distribution $q (z_{i}, ω_{i})$ is a categorical distribution, and $q (z_{i l} = 1, ω_{i c} = 1)$ can be calculated as follows using the term $η_{i c l}$ :

log η_{i c l} = \sum_{j = 1}^{J} \sum_{h = 1}^{2^{\sum_{k} q_{j k}}} g_{j h l} [x_{i j} E (log θ_{j h}) + (1 - x_{i j}) E (log (1 - θ_{j h}))]

+ \sum_{p = 1}^{P} [y_{i p} E (log λ_{p c}) + (1 - y_{i p}) E (log (1 - λ_{p c}))]

+ E (log τ_{c l}) + E (log π_{c}) + const .,

where “ $const .$ ” is a term that does not include parameters. Here, if the random variable W follows a beta distribution with parameters a and b, the expectations of $log W$ and $log (1 - W)$ are

E (log W) = ψ (a) - ψ (a + b),

E (log (1 - W)) = ψ (b) - ψ (a + b),

where $ψ (\cdot)$ is a digamma function. In addition, if a random vector $W = {[W_{1}, \dots W_{m}, \dots W_{M}]}^{⊤}$ follows a Dirichlet distribution with parameter $β = {[β_{1}, \dots, β_{m}, \dots, β_{M}]}^{⊤}$ , the expectation of $log W_{m}$ is

E (log W_{m}) = ψ (β_{m}) - ψ (\sum_{m = 1}^{M} β_{m}) .

Then, $q (z_{i l} = 1, ω_{i c} = 1)$ can be calculated to normalize $η_{i c l}$ as

q (z_{i l} = 1, ω_{i c} = 1) = \frac{η_{i c l}}{\sum_{c = 1}^{C} \sum_{l = 1}^{2^{K}} η_{i c l}} = E (z_{i l} ω_{i c}),

and $q (z_{i l} = 1, ω_{i c} = 1)$ corresponds to the expectation of the product of two latent variables $z_{i l} ω_{i c}$ . This joint probability also leads to marginal probabilities:

q (z_{i l} = 1) = \frac{\sum_{c = 1}^{C} η_{i c l}}{\sum_{c = 1}^{C} \sum_{l = 1}^{2^{K}} η_{i c l}} = E (z_{i l}),

q (ω_{i c} = 1) = \frac{\sum_{l = 1}^{2^{K}} η_{i c l}}{\sum_{c = 1}^{C} \sum_{l = 1}^{2^{K}} η_{i c l}} = E (ω_{i c}) .

Combining these results, the VB algorithm can be summarized as shown in Algorithm 1. This algorithm ensures a monotonic increase in ELBO, and after several iterations, we obtain the optimal variational posteriors.

Algorithm 1: VB Algorithm

Input: Data matrices X and Y; number of latent classesC; Q-matrix; G-matrices; hyper parameters

A^{0}, B^{0}, C^{0}, D^{0}, Γ^{0}, and δ^{0}

; and stopping criterion

ϵ

.
Output: Variational posteriors:

q (z_{i}, ω_{i})

q (θ_{h j})

q (λ_{p c})

q (τ_{c})

, and

q (π)

.
Initialization:

E (z_{i l} ω_{i c}), \forall i, l,

and c, and calculate

E (z_{i l})

and

E (ω_{i c}),

1: Variational-maximization (VM) step: update

q (θ_{j h}), \forall j, h

q (λ_{p c}), \forall p, c

q (τ_{c}), \forall c

, and

q (π)

to calculate Equations 18, 20, 23, and 24.
2: Variational-expectation (VE) step: update

E (z_{i l} ω_{i c}), \forall i, l,

and c and calculate

E (z_{i l})

and

E (ω_{i c})

to calculate Equations 28–30.
3: Calculate ELBO using Equation 14 in Supplementary Material B.
4. If the change in ELBO

< ϵ

, end the algorithm, else return to the VM step.

Note that Algorithm 1 starts from the VM step because initializing

E (z_{i l} ω_{i c})

is more straightforward than setting variational posterior parameters. After the algorithm converges, point estimates, such as the means and variances, are calculated according to the basic results of the beta and Dirichlet distributions. Again, if the random variable W follows a beta distribution with parameters a and b, then

E (W) = a / (a + b)

and

V (W) = a b / {{(a + b)}^{2} (a + b + 1)}

. In addition, if a random vector

W = {[W_{1}, \dots W_{m}, \dots W_{M}]}^{⊤}

follows a Dirichlet distribution with parameter

β

, then

E (W_{m}) = β_{m} / (\sum_{m} β_{m})

and

V (W_{m}) = β_{m} (\sum_{m} β_{m} - β_{m}) / {{(\sum_{m} β_{m})}^{2} (\sum_{m} β_{m} + 1)}

. The estimates of attribute mastery pattern number and latent class of individual i are

l = {argmax}_{l} E (z_{i l})

and

c = {argmax}_{c} E (ω_{i c})

, respectively.

2.3. Gibbs Sampling Algorithm

The results of the VB algorithm can be employed to derive the Gibbs sampling algorithm for two-level DCMs. The full conditional distributions are required to construct the Gibbs sampling algorithm and are easily obtained from the VB algorithm results simply by replacing the expectations with MCMC samples. For example, when the tth iteration of MCMC sample $z_{i l}^{(t)}$ s, whose upper script $(t)$ represents the iteration number, is obtained, the full conditional distribution of $θ_{j h}$ is a beta distribution, and its parameters are

{\begin{matrix} a_{j h}^{(t)} = \sum_{i = 1}^{I} \sum_{l = 1}^{2^{K}} z_{i l}^{(t)} g_{j h l} x_{i j} + a_{j h}^{0}, \\ b_{j h}^{(t)} = \sum_{i = 1}^{I} \sum_{l = 1}^{2^{K}} z_{i l}^{(t)} g_{j h l} (1 - x_{i j}) + b_{j h}^{0} . \end{matrix}

Note that $θ_{j h}$ can be sampled to satisfy the monotonicity constraints using the sampling procedure of Yamaguchi and Templin (2022a). The other full conditional posteriors are listed hereafter.

The full conditional distribution of $λ_{p c}$ is a beta distribution with two parameters (White & Murphy, 2014):

{\begin{matrix} c_{p c}^{(t)} = \sum_{i = 1}^{I} ω_{i c}^{(t)} y_{i p} + c_{p c}^{0}, \\ d_{p c}^{(t)} = \sum_{i = 1}^{I} ω_{i c}^{(t)} (1 - y_{i p}) + d_{p c}^{0} . \end{matrix}

Similarly, the full conditional distributions of $τ_{c}$ and $π$ are Dirichlet distributions with parameters $γ_{c}^{(t)} = {[γ_{c 1}^{(t)}, \dots, γ_{c l}^{(t)}, \dots, γ_{c 2^{K}}^{(t)}]}^{⊤}$ and $δ^{(t)} = {[δ_{1}^{(t)}, \dots, δ_{c}^{(t)}, \dots, δ_{C}^{(t)}]}^{⊤}$ , and the elements of the parameters are

γ_{c l}^{(t)} = \sum_{i = 1}^{I} ω_{i c}^{(t)} z_{i l}^{(t)} + γ_{c l}^{0},

δ_{c}^{(t)} = \sum_{i = 1}^{I} ω_{i c}^{(t)} + δ_{c}^{0} .

Finally, the joint full conditional distribution of z_i and $ω_{i}$ is a categorical distribution with parameter $η_{i c l}^{(t)} / (\sum_{c = 1}^{C} \sum_{l = 1}^{2^{K}} η_{i c l}^{(t)}), \forall c, l$ , where

log η_{i c l}^{(t)} = \sum_{j = 1}^{J} \sum_{h = 1}^{2^{\sum_{k} q_{j k}}} g_{j h l} [x_{i j} log θ_{j h}^{(t)} + (1 - x_{i j}) log (1 - θ_{j h}^{(t)})]

+ \sum_{p = 1}^{P} [y_{i p} log λ_{p c}^{(t)} + (1 - y_{i p}) log (1 - λ_{p c}^{(t)})] + log τ_{c l}^{(t)} + log π_{c}^{(t)} + const .

Assembling these elements, the Gibbs sampling algorithm for two-level DCMs is shown in Algorithm 2.

Algorithm 2: Gibbs Sampling Algorithm

Input: Data matrices X and Y; number of latent classes

C;

Q-matrix; G-matrices; hyper parameters

A^{0}, B^{0}, C^{0}, D^{0}, Γ^{0}, and δ^{0}

; and the number of MCMC iterations T.
Output: MCMC samples of

z_{i}, ω_{i}, θ_{h j}

λ_{p c}

τ_{c}

, and

π

.
Initialization: set

z_{i l}^{(0)}

and

ω_{i c}^{(0)}, \forall i, l,

iteration counter

t = 1

.
1: For all j and h, sample

θ_{j h}^{(t)}

from the beta distribution with parameters

a_{j h}^{(t)}

and

b_{j h}^{(t)}

,which are expressed in Equation 31, satisfying monotonicity constraints.
2: For all p and c, sample

λ_{p c}^{(t)}

from the beta distribution with parameters

c_{p c}^{(t)}

and

d_{p c}^{(t)}

,which are expressed in Equation 32.
3: For all c, sample

τ_{c}^{(t)}

from the Dirichlet distribution with parameter

γ_{c}^{(t)}

,which is expressed in Equation 33.
4: Sample

π^{(t)}

from the Dirichlet distribution with parameter

δ^{(t)}

is expressed in Equation 34.
5:For all

i,

samples

z_{i}^{(t)}

and

ω_{i}^{(t)}

from the categorical distribution with the probability parameters are calculated using Equation 35.
6: If

t = T,

then end the algorithm, else increase the iteration count

t \leftarrow t + 1

and repeat Steps 1–5.

3. Simulation Study

3.1. Simulation Settings

The simulation study aimed to verify the parameter recovery of the developed VB algorithm and Gibbs sampling method. In other words, the objective of the simulation study was to confirm how the proposed estimation algorithms work. The simulation was conducted on a desktop computer with a Windows 10 Pro operating system, Intel Xeon central processing unit with 3.60 GHz and 32.0-GB RAM. The entire simulation code was written in R (Version 4.2.0; R Core Team, 2022). Four factors were manipulated: (1) sample size (200 or 2,000), (2) Q-matrix type (three or four attributes), (3) number of latent classes (three or four), and (4) item quality of X and Y (high or low, which refers to the discrimination power of items for the attribute masteries and latent classes). We also conducted model comparison simulation between an ordinal general DCM and proposed two-level DCM to check whether using additional latent classes improves attribute mastery pattern recovery. This comparative simulation was not the primary research objective of this study; thus, related simulation settings and results are presented in Online Supplementary Material C.

The sample size setting could be small (200) or large (2,000). Sessoms and Henson (2018), in their review of DCM application, noted that sample size is generally large (median = 1,255) with large variations; therefore, we set 2,000 as the large sample size setting. The number of indicators used to define the exogenous latent class was fixed at 24. We controlled the severeness of the condition to change the number of latent classes. However, if we increase the number of indicators according to the increase in the number of latent classes, the information required to distinguish latent classes would increase. Therefore, we must fix the number of indicators of latent classes. Furthermore, the number of latent classes was set as three or four because latent class analysis application studies in an educational setting have sometimes involved four classes (e.g., Lin & Tai, 2015; Toker & Green, 2021). As other examples, Tuominen et al. (2020) analyzed the achievement goal orientation of Finnish sixth- and seventh-grade students and revealed four types of time-invariant classes. Furthermore, Sideridis et al. (2021) applied a multilevel latent class analysis to the achievements of high school students and detected four latent classes connecting GAT science test results and background information about the students. We assumed the same situation as in the previous studies. The three and four latent class situation is also shown in the real data analysis in the next section. The reason we employed 24 indicators is that it is a common multiple of three and four. Twelve is the least common multiple; however, in this case, only three indicators strongly connected to a latent class in the four latent class condition. This might be restrictive; thus, we assumed 24 indicators rather than 12.

Two types of Q-matrices are described in Table 1. The left and right parts of Table 1 show a three-attribute Q-matrix with 20 items and a four-attribute matrix with 22 items, respectively. These numbers of attributes may be small, but the famous Examination for the Certificate of Proficiency in English (e.g., Templin & Hoffman, 2013) data contain three attributes. In another DCM simulation study, Sen and Cohen (2021) also selected three attributes. A situation with four attributes is more challenging than one with three attributes.

Table 1.

Three (Left) and Four (Right) Attributes of Q-Matrices

Item	Attribute			Item	Attribute
Item	1	2	3	Item	1	2	3	4
1	1	0	0	1	1	0	0	0
2	0	1	0	2	0	1	0	0
3	0	0	1	3	0	0	1	0
4	1	0	0	4	0	0	0	1
5	0	1	0	5	1	0	0	0
6	0	0	1	6	0	1	0	0
7	1	0	0	7	0	0	1	0
8	0	1	0	8	0	0	0	1
9	0	0	1	9	1	0	0	0
10	1	1	0	10	0	1	0	0
11	1	0	1	11	0	0	1	0
12	0	1	1	12	0	0	0	1
13	1	1	0	13	1	1	0	0
14	1	0	1	14	1	0	1	0
15	0	1	1	15	1	0	0	1
16	1	1	0	16	0	1	1	0
17	1	0	1	17	0	1	0	1
18	0	1	1	18	0	0	1	1
19	1	1	1	19	1	1	1	0
20	1	1	1	20	1	1	0	1
				21	1	0	1	1
				22	0	1	1	1

The true $Θ$ was generated using the method employed by Yamaguchi and Templin (2022a), in which the correct item response probabilities of all nonmastering item-specific patterns and mastering patterns are represented as g and $1 - s$ . In addition, the number of mastering attributes in the hth item-specific pattern is defined as $#_{j h}$ , and the true correct item response probability is given by $θ_{j h} = g + #_{j h} (1 - s - g) / (\sum_{k} q_{j k})$ . The high item quality conditions were $1 - s = 0.9$ and $g = 0.1$ , and the low item quality conditions were $1 - s = 0.7$ and $g = 0.3$ . In the exogenous latent model measurement part, each lament class had $24 / C$ strong indicators; therefore, the corresponding $λ$ values were generated from $U n i f (s_{l o w e r}, s_{u p p e r})$ . The other $24 (C - 1) / C$ items were weak indicators, and their $λ$ values were obtained from $U n i f (w_{l o w e r}, w_{u p p e r})$ . Based on these considerations, the strong indicator design represented easily interpretable latent classes. By contrast, the weak indicator design represented the low interpretable situation. In other words, the latent class is not clearly defined by the latent classes. The high item quality conditions were set as $s_{l o w e r} = 0.8, and s_{u p p e r} = 0.9$ and $w_{l o w e r} = 0.1, and w_{u p p e r} = 0.2$ , and the low item quality conditions were fixed at $s_{l o w e r} = 0.6, and s_{u p p e r} = 0.8$ and $w_{l o w e r} = 0.2, and w_{u p p e r} = 0.3$ . For example, in the case with $C = 3$ and high quality, the first eight $λ$ values of the first latent class ( $λ_{11}$ to $λ_{81}$ ) were generated from $U n i f (0.8, 0.9)$ , whereas $λ_{91}$ to $λ_{24, 1}$ were obtained from $U n i f (0.1, 0.2)$ .

The true $τ_{c}$ values were obtained from a Dirichlet distribution with a parameter vector of length $2^{K}$ with 1/K as all its values. This choice was made because if the Dirichlet distribution parameters are all 1, then the distribution is uniform; furthermore, smaller parameters of the Dirichlet distribution generate more concentrated probabilities in a few elements than larger parameters. Thus, a latent class connects to a few attribute mastery patterns. However, if the parameters of the Dirichlet distribution are too small, then almost all probabilities are close to 0. The previous trial revealed that using $1 / K$ for all parameters was not too strict but rather appropriately reflected our simulated situation. $π$ was generated from a Dirichlet distribution with a parameter vector of length C, ${[1, 2, \dots, C]}^{⊤}$ . The latent class indicator $Ω$ was generated first, and Z was generated based on $Τ$ and $Ω$ . Using these latent variables and parameters $Θ$ and $Λ$ , the observed X and Y results were generated. Random values were generated for each simulation, and the number of repetitions was 200.

The maximum number of iterations of the VB approach was 1,000, and the convergence criterion was $10^{- 4}$ . The number of iterations of the MCMC algorithm was 4,000, and a single chain was utilized to reduce the risk of label switching. Note that we employed the Geweke (1992) method, which is available in the CODA package (Plummer et al., 2006), and a trace plot for MCMC convergence check in the preliminary study. We judged that 4,000 iterations were sufficient for this simulation. Hyperparameters $A^{0}$ and $B^{0}$ were set to be the same as those in Yamaguchi and Okada (2020a). To avoid label switching, $δ^{0}$ was set to ${[C, C - 1, \dots, 1]}^{⊤}$ . The other hyperparameters were all set to 1. The initial values in the VB method were set as $E (z_{i l} ω_{i c}) = 1 / (C 2^{K}), \forall i, l, c$ . The initial attribute mastery patterns for the Gibbs sampling method were based on the subscore, which is given by $x_{i k}^{s u b s c o} = \sum_{j} x_{i j} q_{j k}$ . The initial $α_{i k}$ was calculated as $α_{i k} = I (x_{i k}^{s u b s c o} > M e d (x_{i k}^{s u b s c o}))$ , where $I (\cdot)$ indicates that its argument is true if its value is 1 and 0 otherwise, and $M e d (x_{i k}^{s u b s c o})$ is the median of the subscore. Finally, the initial $ω_{i c}$ was generated by uniform distortion.

The bias and root-mean-square error (RMSE) were calculated for each parameter. The bias of the mth parameter $ϑ_{m}$ was taken to be $B i a s_{m} = \sum_{t = 1}^{200} ({\hat{ϑ}}_{m t} - ϑ_{m t}^{True}) / 200$ , and its RMSE was obtained from $R M S E_{m} = \sqrt{\sum_{t = 1}^{200} {({\hat{ϑ}}_{m t} - ϑ_{m t}^{True})}^{2} / 200}$ , where ${\hat{ϑ}}_{m t}$ is the estimate of the true parameter $ϑ_{m t}^{True}$ and the subscript t is the replication time. Then, the values were averaged over each parameter set. For example, we calculated the bias or RMSE of $θ_{j h}$ over 200 repetitions and then took the average of this quantity over $Θ$ . The average recovery rates of the elementwise attribute mastery and attribute mastery patterns were also evaluated. For example, the average recovery rate of the kth attribute is $\sum_{t = 1}^{200} \sum_{i = 1}^{I} I ({\hat{α}}_{i k t} = α_{i k t}^{True}) / (I \times 200)$ , where $I (\cdot)$ indicates that its argument is true if it returns 1 and 0 otherwise, ${\hat{α}}_{i k t}$ is the estimated kth attribute of individual i at the tth repetition, and $α_{i k t}^{True}$ is the corresponding true attribute mastery. The true latent class recovery rate was calculated for each class based on the estimated $Ω$ in the same manner. We also calculated the estimation time to check the estimation speed of the developed methods.

3.2. Results

Table 2 summarizes the mean estimation times and standard deviations of the VB and Gibbs sampling algorithms and the ratio of the mean estimated times of both methods. In this simulation, the VB estimation time was affected by the item quality and sample size, where a lower item quality and larger sample size increased the estimation time. For a sample size of 200, the VB estimation was completed within 3 seconds on average. Even with a sample size of 2,000, the estimation took ∼20 seconds at most, even with low item quality. The time taken by the Gibbs sampling algorithm, on the contrary, was mainly affected by the sample size. Specifically, the Gibbs sampling algorithm took approximately 60 seconds and 360–420 seconds when the sample size was 200 and 2,000, respectively. In this simulation study, the VB algorithm was at least ∼20 times faster than the Gibbs sampling method. Under some conditions, such as a sample size of 200 and high item quality, the VB method was more than 100 times faster than the MCMC method. Note that we also checked the mean iteration number of the VB algorithm at the algorithm finishing point for each condition. The result is shown in Online Supplementary Material D. The values were considerably less than the number of maximum iterations (1,000). This indicated that VB estimation was finished before reaching the maximum iteration number. Therefore, we can conclude that VB is considerably faster than the MCMC method.

Table 2.

Estimation Times of Variational Bayesian (VB) and Gibbs Sampling Methods

Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	VB		Gibbs		Gibbs/VB
Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	Mean	(SD)	Mean	(SD)	Gibbs/VB
200	Three attributes	Three	High	0.305	(0.105)	51.245	(2.374)	168.009
		Three	Low	1.398	(0.660)	51.220	(2.134)	36.631
		Four	High	0.344	(0.162)	52.871	(2.393)	153.769
		Four	Low	1.520	(0.715)	52.549	(2.180)	34.574
	Four attributes	Three	High	0.401	(0.108)	58.448	(2.346)	145.891
		Three	Low	2.346	(1.296)	58.082	(2.317)	24.756
		Four	High	0.494	(0.181)	60.945	(2.606)	123.249
		Four	Low	2.728	(1.519)	60.906	(2.441)	22.326
2,000	Three attributes	Three	High	3.145	(3.302)	367.194	(13.150)	116.772
		Three	Low	11.171	(4.529)	365.144	(12.457)	32.687
		Four	High	3.780	(2.531)	376.420	(16.839)	99.587
		Four	Low	11.422	(5.053)	373.769	(21.694)	32.723
	Four attributes	Three	High	5.700	(4.947)	399.276	(19.673)	70.048
		Three	Low	20.768	(6.417)	397.060	(43.274)	19.119
		Four	High	5.738	(3.628)	417.557	(19.498)	72.775
		Four	Low	20.530	(4.732)	419.173	(22.654)	20.417

Table 3 summarizes the biases and RMSEs of correct item response probability parameter Θ. The VB and Gibbs sampling methods exhibited similar small biases, where a sample size of 2,000 produced smaller biases than that of 200. The sample size of 200 and lower item quality yielded small negative biases (–0.011 to –0.018). However, the other cases exhibited biases of less than $- 0.005$ and sufficiently small parameter recovery. The RMSEs of the VB method with a sample size of 200 were larger than those of the Gibbs sampling method, but the relationship was the opposite in the 2,000 sample size case. However, the RMSEs were generally smaller than 0.05, indicating that the parameter estimates of the two Bayesian methods were stable.

Table 3.

Biases and Root-Mean-Square Errors (RMSEs) of Correct Item Response Probability Parameter Θ of Variational Bayesian (VB) and Gibbs Sampling Methods

Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	VB		Gibbs
Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	Bias	(RMSE)	Bias	(RMSE)
200	Three attributes	Three	High	−.005	(.011)	−.005	(.006)
		Three	Low	−.016	(.026)	−.018	(.008)
		Four	High	−.004	(.011)	−.004	(.006)
		Four	Low	−.015	(.026)	−.015	(.007)
	Four attributes	Three	High	−.003	(.011)	−.003	(.005)
		Three	Low	−.012	(.027)	−.012	(.008)
		Four	High	−.004	(.010)	−.004	(.005)
		Four	Low	−.011	(.027)	−.011	(.008)
2,000	Three attributes	Three	High	−.001	(.003)	−.002	(.066)
		Three	Low	−.002	(.009)	−.003	(.009)
		Four	High	.000	(.002)	−.003	(.051)
		Four	Low	−.003	(.006)	−.002	(.009)
	Four attributes	Three	High	.000	(.003)	−.001	(.041)
		Three	Low	−.002	(.009)	−.001	(.012)
		Four	High	.000	(.002)	−.001	(.032)
		Four	Low	−.002	(.007)	−.002	(.008)

Table 4 presents the biases and RMSEs of relationship parameter $Τ$ . Both the VB and Gibbs sampling methods provide unbiased estimates of $Τ$ in this simulation. Low-quality items exhibited larger RMSEs than high-quality items, but the RMSEs remained less than 0.04 and were sufficiently small even when the sample size was 200. The RMSEs with a sample size of 2,000 were less than those with a sample size of 200.

Table 4.

Biases and Root-Mean-Square Errors (RMSEs) of Relationship Parameter $Τ$ of Variational Bayesian (VB) and Gibbs Sampling Methods

Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	VB		Gibbs
Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	Bias	(RMSE)	Bias	(RMSE)
200	Three attributes	Three	High	.000	(.016)	.000	(.012)
		Three	Low	.000	(.031)	.000	(.026)
		Four	High	.000	(.017)	.000	(.013)
		Four	Low	.000	(.039)	.000	(.032)
	Four attributes	Three	High	.000	(.008)	.000	(.007)
		Three	Low	.000	(.020)	.000	(.018)
		Four	High	.000	(.012)	.000	(.010)
		Four	Low	.000	(.025)	.000	(.023)
2,000	Three attributes	Three	High	.000	(.009)	.000	(.026)
		Three	Low	.000	(.006)	.000	(.011)
		Four	High	.000	(.008)	.000	(.021)
		Four	Low	.000	(.011)	.000	(.015)
	Four attributes	Three	High	.000	(.005)	.000	(.005)
		Three	Low	.000	(.005)	.000	(.007)
		Four	High	.000	(.004)	.000	(.005)
		Four	Low	.000	(.006)	.000	(.007)

Table 5 displays the biases and RMSEs of population mixing parameter $π$ . The RMSEs were almost the same between the VB and Gibbs sampling methods. Most biases were less than 0.01, but slightly larger biases of 0.014–0.018 were evident in the cases of 200 sample size and four latent classes. The biases obtained with the larger sample size conditions were satisfactorily small. The VB method tended to have larger biases than the Gibbs sampling method. Table 6 provides the biases and RMSEs of parameter $Λ$ . As shown in the table, a larger sample size reduced the biases in the four latent cases. The RMSEs of the VB method were larger than those of the Gibbs sampling method. Thus, the former was less stable than the latter, but most RMSEs were less than 0.04.

Table 5.

Biases and Root-Mean-Square Errors (RMSEs) of Population Mixing Parameter of Variational Bayesian (VB) and Gibbs Sampling Methods

Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	VB		Gibbs
Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	Bias	(RMSE)	Bias	(RMSE)
200	Three attributes	Three	High	.007	(.024)	.002	(.023)
		Three	Low	.003	(.020)	.001	(.018)
		Four	High	.008	(.008)	.001	(.005)
		Four	Low	.010	(.016)	.006	(.017)
	Four attributes	Three	High	.008	(.009)	.000	(.001)
		Three	Low	.001	(.003)	.001	(.008)
		Four	High	.011	(.016)	.001	(.006)
		Four	Low	.013	(.018)	.004	(.013)
2,000	Three attributes	Three	High	.003	(.005)	.000	(.000)
		Three	Low	.000	(.000)	.000	(.002)
		Four	High	.002	(.003)	.000	(.000)
		Four	Low	.002	(.005)	.002	(.005)
	Four attributes	Three	High	.005	(.005)	.000	(.000)
		Three	Low	.000	(.000)	.000	(.000)
		Four	High	.005	(.005)	.000	(.000)
		Four	Low	.000	(.000)	.000	(.000)

Table 6.

Biases and Root-Mean-Square Errors (RMSEs) of Response Probability Parameter Λ for Variational Bayesian (VB) and Gibbs Sampling Methods

Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	VB		Gibbs
Sample Size	Q-Matrix	Number of Latent Classes	Item Quality	Bias	(RMSE)	Bias	(RMSE)
200	Three attributes	Three	High	.009	(.032)	.008	(.019)
		Three	Low	.008	(.023)	.007	(.018)
		Four	High	.016	(.039)	.016	(.021)
		Four	Low	.014	(.035)	.014	(.024)
	Four attributes	Three	High	.009	(.037)	.008	(.015)
		Three	Low	.009	(.024)	.008	(.016)
		Four	High	.018	(.043)	.014	(.021)
		Four	Low	.017	(.037)	.015	(.024)
2,000	Three attributes	Three	High	.003	(.021)	.001	(.003)
		Three	Low	.002	(.011)	.001	(.007)
		Four	High	.004	(.018)	.003	(.004)
		Four	Low	.004	(.016)	.002	(.008)
	Four attributes	Three	High	.003	(.026)	.001	(.004)
		Three	Low	.002	(.010)	.001	(.007)
		Four	High	.004	(.020)	.002	(.005)
		Four	Low	.005	(.013)	.003	(.006)

Figure 2 depicts the recovery of attribute mastery when the Q-matrix has four attributes. The attribute and attribute mastery patterns are generally well recovered, and the two estimation methods provide almost identical results. The lower quality items reduce the attribute recovery rate. Figure 3 presents the latent class recovery results obtained when there are four latent classes. Again, the latent class is well recovered with both methods. Note that the results are almost the same when the Q-matrix has three or four attributes, as well as when there are three or four latent classes. Therefore, the attribute recovery results in the case of a Q-matrix with three attributes, and the latent class recovery results when three latent classes are omitted.

Figure 2.

Attribute mastery recovery rate in the case of four attributes.

Figure 3.

Latent class recovery rate in the case of four classes.

In summary, the VB method is faster than the Gibbs sampling method, and the two estimation methods are unbiased. The Gibbs sampling method yields more precise and stable results than the VB method. However, both methods provide satisfactory recovery and effectively recover the individual latent variables that are attributes and latent classes. Nevertheless, both methods have advantages and disadvantages; thus, the estimation method should be selected considering these aspects.

4. Real Data Analysis

4.1. Data Analysis Settings

Fourth-grade data from TIMSS 2007 were employed to demonstrate the proposed two-level DCM. Specifically, 25 items from Booklets 4 and 5 of the TIMSS 2007 fourth-grade mathematics assessment were selected for the DCM measurement, which is also consistent with the previous studies (Lee et al., 2011; Yamaguchi & Okada, 2018). The Q-matrix included in the DCM package is summarized in Table S4 in Online Supplementary Material D (George et al., 2016). Three domains of test items were used as attributes for simplicity: number, geometric shapes and measures, and data display. The original set defined by Lee et al. (2011) contained 15 attributes, but it was too complex; therefore, we reduced the number of attributes. We selected country data from the United States and Canada, including Massachusetts, Minnesota, Alberta, British Columbia, Ontario, and Quebec.

As the latent classes, we selected 19 items from the student questionnaires on mathematics in school obtained from TIMMS 2007 Supplemental 1 pp. 16–18 (Foy & Olson, 2009). The items were from AS4MAWEL (“I usually do well in mathematics”) to AS4MHCOM (“I use a computer”). The first eight items were assessed with a four-point Likert-type scale (1 = agree a lot, 2 = agree a little, 3 = disagree a little, and 4 = disagree a lot), and the remaining 12 items were analyzed on another four-point Likert-type scale with different labels (1 = every or almost every lesson, 2 = about half the lessons, 3 = some lessons, and 4 = never). We revised the original four-point scale using dichotomous variables to treat 1 and 2 as 1, and 3 and 4 as 0, respectively. Therefore, a recorded value of 1 means agreement regarding attitude toward mathematics or performing mathematical activities for more than half of the lesson. Note that dichotomizing Likert items leads to loss of information. We do not believe the procedure is always appropriate. Without dichotomization, however, the interpretation of latent classes based on the estimated item response probabilities became complex and the graph became very messy. Therefore, for the sake of simplicity, we dichotomized the questionnaire items in this study. However, the latent class with the polytomous response could be incorporated into the two-level DCMs. The individuals with missing values were eliminated, and the total sample size was 1,061.

Preliminary latent class analysis was conducted to determine the number of latent classes using the poLCA package (Linzer & Lewis, 2011). In addition, we fit two to six classes in the two-level DCMs to determine the number of latent classes. The maximum number of iterations in the VB and MCMC approaches was 500 and 8,000, respectively. The other estimation settings, such as the hyperparameter settings, were the same as those in the simulation study. Both estimation methods were compared in terms of whether they provided similar results. Furthermore, we determined the latent class results and relationships between the latent classes and attribute mastery patterns. Data analysis codes and data are available from the Open Science Framework (OSF: https://osf.io/6nxtc/) page.

4.2. Results

Table 7 presents the fit indices of the preliminary latent class analysis results. Bayesian information criterion (BIC) indicates that four classes were the best, and the other indices, such as Akaike information criterion (AIC) or G ², continued to decrease with an increasing number of latent classes. However, the changes in these values became small after four classes.

Table 7.

Preliminary Analysis to Determine the Number of Latent Classes

Number of Latent Classes	Log-Likelihood	Number of Parameters	Akaike Information Criterion	Bayesian Information Criterion	G ²
2	−10,182.73	39	20,443.45	20,637.16	6,799.20
3	−9,874.63	59	19,867.26	20,160.31	6,183.01
4	−9,757.86	79	19,673.71	20,066.10	5,949.46
5	−9,700.92	99	19,599.85	20,091.58	5,835.60
6	−9,652.10	119	19,542.20	20,133.27	5,737.95

We also calculated AIC and BIC of the two-level DCMs based on the VB estimation results, as shown in Table 8. AIC indicated that having four classes was the best, while BIC indicated three. ELBO indicated a similar tendency in the log-likelihood shown in Table 7; the change in ELBO became small after four classes. Therefore, combining the results from Tables 7 and 8, three or four classes are possible. We also considered the interpretability of the latent class and thus selected four classes in this study. The possibility of different numbers of classes will be discussed later.

Table 8.

Information Criteria of the Two-Level Diagnosis Classification Models Based on the Variational Bayesian Results

Number of Classes	Evidence Lower Bound	Number of Parameters	Akaike Information Criterion	Bayesian Information Criterion
2	−24,925.62	119	50,089.23	50,680.30
3	−24,679.41	146	49,650.81	50,375.99
4	−24,587.51	173	49,521.03	50,380.31
5	−24,569.77	200	49,539.53	50,532.93
6	−24,560.60	227	49,575.20	50,702.70

Next, we specified the four latent classes in two-level DCMs and estimated model parameters using the VB and Gibbs sampling methods. The trace plot obtained after a burn-in period on OSF indicated that no systematical trend existed, and we judged that the MCMC chains had converged. Figure 4 provides a collection of scatter plots of the parameter estimates obtained with VB and Gibbs sampling methods for four parameter sets. All correlations between the two methods are greater than .99; thus, the parameter estimates resulting from both methods are almost the same. Therefore, the parameter estimation results of the VB method are reported hereafter.

Figure 4.

Scatter plots of variational Bayesian and Gibbs sampling methods (Markov chain Monte Carlo) for Θ (upper left panel), Λ (upper right panel), Τ (lower left panel), and π (lower right panel).

Figure 5 presents the estimates of latent class response probability parameter $Λ$ in the four classes case. The first class shows high response probabilities for the 1st (AS4MAWEL: “I usually do well in mathematics”), 7th (AS4MABOR: “Mathematics is boring”), 14th (AS4MHMWP: “I memorize how to work problems”), and 17th (AS4MHWPO: “I work problems on my own”) items especially. However, the first class shows low response probabilities for the 2nd (AS4MAMOR: “I would like to do more mathematics in school”), 4th (AS4MAENJ: “I enjoy learning mathematics”), 8th (AS4MALIK: “I like mathematics”), and 11th (AS4MHMCL:“I measure things in the classroom and around the school”) items. These response probabilities indicate that the class feels efficiency for mathematics but is boring and disliked.

Figure 5.

Estimated parameter Λ of four latent classes with variational Bayesian estimation. Note. Each number represents the corresponding class.

The second class exhibits moderate response probabilities for almost all items. The second class produces relatively high probabilities for the 3rd (AS4MACLM: “Mathematics is harder for me than for many of my classmates”) and 5th (AS4MANOT: “I am just not good at mathematics”) items and lower for the 1st, 6th (AS4MAQKY: “I learn things quickly in mathematics”), 9th (AS4MHASM: “I practice adding, subtracting, multiplying, and dividing without using a calculator”), 14th, and 17th items. This class hated mathematics and did not actively work in the classroom.

The third and fourth classes tend to have high probabilities for the first, second, fourth, sixth, and eighth items, implying that both groups like mathematics. However, these classes show discrepancies in the items after the ninth item. The third class tends to have lower probabilities for the latter half of the items than the fourth class. Large discrepancies are evident for items such as the 10th (AS4MHWFD: “I work on fractions and decimals”), 11th (AS4MHMCL: “I measure thing in the classroom and around the school”), and 15th (AS4MHWSG: “I work with other students in small groups”) items. These items indicate engagement with the mathematical activities in the class. In summary, the fourth class likes mathematics and engages in mathematical activities, and the third class also likes mathematics but does not engage in mathematical activities.

Figure 6 depicts the estimates of relationship matrix $Τ$ . The attribute mastery patterns, such as only mastering the number attribute (100), only mastering geometric shapes and measures (010), and mastering number and Data and display (101), were not associated with the estimated latent classes. The first class, which had self-efficacy for math but finds that mathematics is boring, tends to master number and geometric shapes and measures (110), only master data and display (001), or master all three attributes (111). In addition, pattern (011), which refers to mastering geometric shapes and measures and data and display, and pattern (001) also showed similar value as pattern (110), that is, fourth. This result indicates that the first class comprised two types of students: those who had high mathematical abilities and found the class too easy and were bored by it, and those who had illusional efficacy for mathematics but did not have the actual skills. The second class, which did not like mathematics, had strong connections to not mastering any of the attributes (000) and only mastering data and display (001). These findings were understandable because the class members dislike mathematics and can only master a few attributes.

Figure 6.

Estimates of relationship matrix Τ for four latent classes with variational Bayesian estimation. Note. Each number represents the corresponding class.

The third class liked mathematics but did not engage in mathematical activities and was strongly connected to mastering all three attributes (111) and only mastering data and display (001). These classes mastered at least one attribute, and many of them master all attributes. The fourth class liked mathematics and engaged in mathematical activities and tended to not master any of the attributes (000), only master data and display (001), or master all three attributes (111). This finding was interesting because mathematical activities in the classroom might not have a strong effect on mastering attributes. Another interpretation was that the mathematical activities in the classroom might be appropriate for students with at most moderate mathematical abilities, and the students who had to engage in the activities did not have high mathematical skills.

In summary, students who felt efficacy or liked mathematics tended to master more attributes, and engagement in mathematics activities might not affect attribute mastering. Interestingly, liking math and being actively engaged in classroom activities were associated with opposite attribute mastery statuses. Subjective engagement might not affect attribute mastery. Furthermore, students who disliked mathematics tended to have less attribute mastery, and it was effective to intervene in this class.

5. Discussion and Conclusion

Two-level DCMs and VB and Gibbs sampling methods were developed. A simulation study showed that both estimation methods provide sufficiently accurate parameter recovery, but the Gibbs sampling method is slightly better than the VB method. However, the VB method is much faster than the Gibbs sampling method. A real data analysis example involving the application of these methods to TIMMS 2007 fourth-grade mathematics data demonstrated how external latent classes are related to attribute mastery patterns. The connection between the external latent class and attribute mastery pattern could reveal the attitudes of students toward mathematics, and engagement in mathematical activities was connected to different attribute mastery patterns. This information can be employed to make teaching plans that consider not only attribute mastery patterns but also factors such as mathematical attitudes.

The simulation study demonstrated that both the proposed estimation methods correctly recovered the model parameters. However, the Gibbs sampling method provided smaller biases and RMSEs than the VB method under some conditions, although the estimation speed of the VB method was much faster than that of Gibbs sampling. Based on these results, the VB method is appropriate for model exploration, and the Gibbs sampling method should be employed in the final estimation step. The estimation methods can be flexibly changed based on the purpose of data analysis and the available computational environment.

The proposed estimation methods were both Bayesian estimation methods. Maximum likelihood estimation is also an alternative estimation method. The EM algorithm can be derived for the ML method because MAP estimation can also be easily derived from the complete data likelihood shown in Equation 12, as discussed in Online Appendix A. Furthermore, the ML method is available via general latent variable modeling software such as Mplus (Muthén & Muthén, 1998–2017) or Latent Gold (Vermunt & Magidson, 2005). Using general latent variable modeling software also extends two-level DCMs to include external covariates. These covariates can be employed to explain the relationship probabilities between latent classes and attribute mastery patterns. One limitation of the proposed Bayesian estimation is that it does not include additional external covariates. Hence, more flexible estimation methods are required in future research.

A two-level DCM is related to various models and can easily be extended. If the measured part of the DCM is the usual latent class model, then two-level DCMs become two-level latent class models (Miyazaki et al., 2007). Furthermore, if the exogenous latent classes are DCMs in the cross-sectional data collection case, then the model represents the strength of the connection between two separate attribute mastery patterns measured by two diagnostic assessments. This model can be used to explore two sets of attribute mastery patterns. For example, the relationships between mathematical attributes and English reading skill sets can be explored simultaneously. Note that there is an arbitrariness as to which attribute set is exogenous when there are no strong theoretical assumptions.

As an extension of these two cross-sectional DCMs, the latent transition or hidden Markov types of longitudinal DMCs (e.g., Chen et al., 2017; Madison & Bradshaw, 2018; Wang, 2021; Yamaguchi & Martinez, in press) can be introduced. Longitudinal DCMs include both endogenous and exogenous latent classes defined by DCM measures in repeated measurement situations. If the diagnostic measurements are performed more than twice, these models can address long-term changes in individual attribute mastery status.

Note that in this example, the exogenous latent class was defined by binary items; however, this measurement model can be changed. For example, polytomous category items or continuous indicators can be employed. Furthermore, general mixture models (e.g., McLachlan & Peel, 2000) can also be employed to define latent classes, which are also called mixture components in the mixture model context. In addition, external continuous observed/latent variables can be included to represent the strength of the connection between the latent class and attribute mastery patterns.

We need several notes on the parameterization of the two-level DCMs. The size of the $Τ$ matrix becomes large if the number of attributes $K$ is medium to large. In modeling each attribute mastery probability, we assume conditional independence among attributes, just as in conventional explanatory DCMs. If the assumption is appropriate, the two-level DCMs will reduce to a more parsimonious form that is a special case of explanatory DCMs. The attribute patternwise model includes such a model because current two-level DCMs do not need the conditional independence assumption. Furthermore, marginalization of the $Τ$ matrix will produce each attribute level connection. In other words, attribute mastery pattern results can be converted into the mastery status of each attribute. In this regard, this patternwise model yields considerable information, thereby having an edge over previous explanatory DCMs. The model formulation can be over parameterization, but there are prior distributions on the $τ_{c}$ vectors. The priors have regularization effects; thus, the Bayesian formulation helps avoid overfitting in case the maximum likelihood estimation does not work. In summary, the attribute patternwise model employed in this study has richer information than the attribute elementwise modeling.

In addition, we need to mention the model specification to understand the structural parameter $τ$ . We can use diagnostic information to infer latent classes. In our model formulation, latent classes are conditioning on attribute mastery probability. However, inversely, we can assume that the attribute mastery patterns are explanatory variables on the latent classes. If a researcher is interested in how attribute mastery patterns affect the latent classes, we can use the latter relationship. Assume the case that students take a diagnostic assessment test at the start of a semester and they also take a mental health survey for adaptation by their school. In such a case, a specific attribute mastery pattern may be able to predict a mental health or adaptation pattern and teachers may use such findings to prevent student drop outs.

As an additional note, if the purpose is only to diagnose students’ knowledge status, the external latent classes may not be necessary. It is important to carefully analyze the meaning of “necessary” here. If we want to understand the relationship between the attribute mastery patterns and outside elements of a diagnostic assessment, the latent classes may be informative. The external information sometimes helps improve the estimation precision of attribute mastery patterns, which was confirmed using additional simulation, but this point is not the primary research topic in this study. We introduced the latent classes to surmise and deeply understand the attribute mastery patterns via the relationship probability parameter τs. The relationship probabilities help understand which latent class is strongly connected to which attribute mastery pattern. Such information is not provided by ordinal DCMs framework; thus, we believe that the latent classes are necessary.

In addition, if the attribute mastery patterns determine the latent classes, the model specification in this study is wrong. In such a case, the relationship matrix is biased because the current model assumes $\sum_{l = 1}^{2^{K}} τ_{l c} = 1$ , but the misspecified model assumes $\sum_{c = 1}^{C} τ_{l c} = 1$ . These different constraints can cause a bias in the structural part. Therefore, it is important to carefully choose which latent variables explain the other in the two-level DCMs.

The interpretation of the $τ$ parameter needs to be made carefully. In the real data analysis, the $τ$ parameters were interpretated based on relative comparison. If a data analyst has a strong opinion about the meaningful cut-off value on the $τ$ parameters, they can employ such values. However, in our data analysis setting, there was no strong opinion on the cut-off values. Therefore, we descriptively interpretated the $τ$ parameters instead of setting strict cut-off values. One option is to set the values to $1 / 2^{K},$ and if $τ$ is greater (or less) than the value, we may conclude the relationship as being stronger (or weaker) than the chance level. However, attributer mastery patterns themselves are not uniformly distributed; thus, we need to check the appropriateness of this cut-off criterion. Furthermore, the cut-off values depend on the context of the test use and careful consideration is required.

The identification of the two-level DCMs must also be discussed. Swapping of attribute mastery and switching of latent class labels may occur if identification conditions are not satisfied. For identification of diagnostic model parts, we need monotonicity constraints on the correct item response probabilities (e.g., Henson et al., 2009; Yamaguchi & Templin, 2022a). In our model formulation, monotonicity constraints were not employed. Instead, our MCMC method satisfies monotonicity constraints for each iteration. However, the VB method does not strictly satisfy monotonicity constraints, but the prior means were set to satisfy monotonicity relationships. Furthermore, one method to prevent label switching of the latent class part is to order parameter constraints on $λ$ s. Instead, to avoid label switching, $δ^{0}$ was set to ${[C, C - 1, \dots, 1]}^{⊤}$ ; however, this cannot completely prevent label switching, especially in each iteration of an MCMC method without constraints. However, our simulation provided satisfactory recovered parameters; thus, we believe that these methods empirically work. In future, more detailed theoretical properties of prior distribution for identification are required.

In addition to the above discussion, one important aspect of model parameter identifiability is the completeness of the Q-matrix: The identity matrix whose size is $K \times K$ should be contained (Chiu et al., 2009). This is a necessary condition for model identification in cases such as the DINA model (Xu & Zhang, 2016). In addition to the completeness of Q-matrix, Xu (2018) established sufficient condition of model parameter for restricted latent class family, which is a wider class of DCMs. The two conditions are: (1) the Q-matrix has two identity matrix parts after row permutation and (2) the remaining part of the Q-matrix should include an identity matrix part. In our real data analysis, the Q-matrix employed was not complete, so we must interpret the results with caution.

The identifiability of DCM parameters is a very important topic not only from a theoretical perspective but also in real data analysis, and the topic has been actively studied for the past decade. Related work by Xu and Shang (2017) proved identifiability conditions for both model parameters and Q-matrix elements. Identifiability conditions in the context of Q-matrix estimation is one important research topic within the research field of DCMs. Moreover, Gu and Xu (2020) extended the identifiability conditions proposed by Xu (2017) under correct model specification and provided strict and partial identifiability conditions that are relatively easily checked. In addition, Chen et al. (2020) modified the model formulation of DCMs to directly utilize the Q-matrix and discussed generic identifiability conditions that are milder than strict identifiability. Culpepper (2023) further relaxed the identifiability condition that requires at least two identity sub-Q-matrices. These studies primarily focused on cross-sectional DCMs. Liu et al. (2023) established strict and genetic identifiability conditions for hidden Markov type DCMs, including the estimation of Q-matrix (see Theorem 3); thus, this work may provide insights into the identifiability of two-level DCMs. Although careful consideration is necessary, based on the Liu et al. (2023), the two-level DCMs model parameters may be identifiable if the population mixing parameter $π$ is positive, the relationship matrix $Τ$ is full rank, the Q-matrix is complete, and the test includes at least one item with a $q$ -vector that generates different item response vector probabilities for an arbitral pair of attribute mastery patterns.

In application, two-level DCMs offer different perspectives on educational intervention with ordinal DCMs. Ordinal DCMs have been employed to improve individual learning using diagnostic feedback based on attribute mastery probabilities or patterns, which is a type of educational intervention. Furthermore, ordinal DCMs can be used not only for individual self-learning but also class room teaching situation. For both situations, the attribute mastery probabilities and patterns are useful information to understand an individual learner’s cognitive weaknesses and strengths. However, diagnostic information may not be sufficient in a class room setting, and it does not tell the reason why students do not master the assumed attributes. Furthermore, it may be naive to think that an educational intervention is equally effective for students who belong to the same attribute mastery pattern. This is an extensive consideration of well-known aptitude–treatment interaction (e.g., Cronbach & Webb, 1975; Snow, 1991). In other words, the effect of intervention can be differed in a specific attribute mastery pattern depending on the variability of their background situations.

The two-level DCMs can provide information about attribute mastery status and other learning activity and individual traits. From the real data analysis, two-level DCMs revealed that students who have already been actively engaged in math class activities do not always master the attributes. This means that some students have already joined the math class but face several difficulties in mastering some attributes. For such students, encouraging them to join the class room activities may not be effective for improving their mathematical abilities. Furthermore, the real data analysis also indicated that some students dislike mathematics, do not engage in the activity, and do not master the attributes. For such students, it may be a good choice to encourage them to join activities and let them enjoy mathematics. These latent classes indicated the variety in non-mastering status of attributes. In a classroom setting, it is difficult to customize the lecture for each student, and there is a need to focus on latent classes that may be most effectively improved. In summary, two-level DCMs can be employed to explore the variability in an attribute mastery status using external information.

Another future research topic relates to determining the latent classes in the two-level DCMs. In the empirical data analysis, we fitted both classical latent class models and proposed two-level DCMs according to various numbers of latent classes. The results indicated that the diagnosis item might change the number of classes. However, this may be the case for only this dataset. Therefore, we need to confirm how the number of latent classes is affected by considering the diagnosis items in a future simulation research.

Supplemental Material

Supplemental Material, sj-docx-1-jeb-10.3102_10769986231173594 - Bayesian Analysis Methods for Two-Level Diagnosis Classification Models

Supplemental Material, sj-docx-1-jeb-10.3102_10769986231173594 for Bayesian Analysis Methods for Two-Level Diagnosis Classification Models by Kazuhiro Yamaguchi in Journal of Educational and Behavioral Statistics

Footnotes

Author’s Note

The data analysis code is available in the Open Science Framework page: .

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 work was supported by Japan Society for the Promotion of Science KAKENHI 19H00616, 20H01720, 21H00936, and 22K13810.

ORCID iD

Kazuhiro Yamaguchi

References

Bishop

(2006). Pattern recognition and machine learning. In Pattern recognition. Springer. https://doi.org/10.1641/B580519

Blei

D. M.

Kucukelbir

McAuliffe

J. D.

(2017). Variational inference: A review for statisticians. Journal of the American Statistical Association, 112, 859–877. https://doi.org/10.1080/01621459.2017.1285773

Brooks

Gelman

Jones

G. L.

Meng

X.-L.

(2011). Handbook of Markov chain Monte Carlo. CRC Press.

Chen

de la Torre

(2014). A procedure for diagnostically modeling extant large-scale assessment data: The case of the programme for international student assessment in reading. Psychology, 5, 1967–1978. https://doi.org/10.4236/psych.2014.518200

Chen

Culpepper

S. A.

Wang

Douglas

(2017). A hidden Markov model for learning trajectories in cognitive diagnosis with application to spatial rotation skills. Applied Psychological Measurement, 42, 5–23. https://doi.org/10.3102/0146621617721250

Chen

Culpepper

Liang

(2020). A sparse latent class model for cognitive diagnosis. Psychometrika, 85(1), 121–153. https://doi.org/10.1007/s11336-019-09693-2

Chiu

C. Y.

Douglas

J. A.

(2009). Cluster analysis for cognitive diagnosis: Theory and applications. Psychometrika, 74(4), 633–665. https://doi.org/10.1007/s11336-009-9125-0

Collins, L. M. & Lanza, S. T. (2009). Latent class and latent transition analysis: With applications in the social, behavioral and health sciences. Wiley.

Cronbach

L. J.

Webb

(1975). Between-class and within-class effects in a reported aptitude treatment interaction: Reanalysis of a study by G. L. Anderson. Journal of Educational Psychology, 67(6), 717–724. https://doi.org/10.1037/0022-0663.67.6.717

10.

Culpepper

S. A.

(2023). A note on weaker conditions for identifying restricted latent class models for binary responses. Psychometrika, 88, 158–174. https://doi.org/10.1007/s11336-022-09875-5

11.

De Boeck

(2004). Explanatory item response models: A generalized linear and nonlinear approach. Springer Science & Business Media.

12.

de la Torre

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

13.

Foy

Olson

J. F.

(2009). TIMSS2007 user guide for the international database; Supplement 1. International version of the TIMSS 2007 background and curriculum questionnaires. TIMSS and PIRLS International Study Center.

14.

George

A. C.

Robitzsch

Kiefer

Groß

Ünlü

(2016). The R package CDM for cognitive diagnosis models. Journal of Statistical Software, 74, 1–24. https://doi.org/10.18637/jss.v074.i02

15.

Geweke

(1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bernardo

J. M.

Berger

J. O.

Dawid

A. P.

Smith

A. F. M.

(Eds.), Bayesian Statistics (4th ed., pp. 169–193). Oxford University Press.

16.

(2020). Partial identifiability of restricted latent class models. The Annals of Statistics, 48(4), 2082–2107. https://doi.org/10.1214/19-AOS1878

17.

Hagenaars

J. A.

McCutcheon

A. L.

(Eds.). (2002). Applied latent class analysis. Cambridge University Press. https://doi.org/10.1017/CBO9780511499531

18.

Henson

R. A.

Templin

J. L.

Willse

J. T.

(2009). Defining a family of cognitive diagnosis models using log-linear models with latent variables. Psychometrika, 74, 191–210. https://doi.org/10.1007/S11336-008

19.

Jeon

Rijmen

Rabe-Hesketh

(2017). A variational maximization–maximization algorithm for generalized linear mixed models with crossed random effects. Psychometrika, 82, 693–716. https://doi.org/10.1007/s11336-017-9555-z

20.

Lazarsfeld

P. F.

Henry

N. W.

(1969). Latent structure analysis. Houghton Mifflin.

21.

Lee

Y.-S.

Park

Y. S.

Taylan

. (2011). A cognitive diagnostic modeling of attribute mastery in Massachusetts, Minnesota, and the U.S. national sample using the TIMSS2007. International Journal of Testing, 11, 144–177. https://doi.org/10.1080/15305058.2010.534571

22.

Lin

S. W.

Tai

W. C.

(2015). Latent class analysis of students’ mathematics learning strategies and the relationship between learning strategy and mathematical literacy. Universal Journal of Educational Research, 3(6), 390–395. https://doi.org/10.13189/ujer.2015.030606

23.

Linzer

D. A.

Lewis

J. B.

(2011). poLCA: An R package for polytomous variable latent class analysis. Journal of Statistical Software, 42(10), 1–29. https://doi.org/10.18637/jss.v042.i10

24.

Liu

Culpepper

S. A.

Chen

(2023). Identifiability of hidden Markov models for learning trajectories in cognitive diagnosis. Psychometrika. https://doi.org/10.1007/s11336-023-09904-x

25.

Madison

M. J.

Bradshaw

L. P.

(2018). Assessing growth in a diagnostic classification model framework. Psychometrika, 83, 963–990. https://doi.org/10.1007/s11336-018-9638-5

26.

McLachlan

Peel

(2000). Finite mixture models. John Wiley & Sons.

27.

Miyazaki

Tamura

Shigemasu

(2007). Two-level latent class analysis with Bayesian inference. Kodo Keiryogaku (The Japanese Journal of Behaviormetrics), 34, 21–31. https://doi.org/10.2333/jbhmk.34.21

28.

Nakajima

Watanabe

Sugiyama

(2019). Variational Bayesian learning theory. Cambridge University Press. https://doi.org/10.1017/9781139879354

29.

Ormerod

J. T.

Wand

M. P.

(2010). Explaining variational approximations. The American Statistician, 64, 140–153. https://doi.org/10.1198/tast.2010.09058

30.

Park

Y. S.

Lee

Y.-S. S.

(2014). An extension of the DINA model using covariates: Examining factors affecting response probability and latent classification. Applied Psychological Measurement, 38, 376–390. https://doi.org/10.3102/0146621614523830

31.

Park

Y. S.

Xing

Lee

(2018). Explanatory cognitive diagnostic models: Incorporating latent and observed predictors. Applied Psychological Measurement, 42, 376–392. https://doi.org/10.3102/0146621617738012

32.

Plummer

Best

Cowles

Vines

(2006). CODA: Convergence diagnosis and output analysis for MCMC. R News, 6, 7–11. https://journal.r-project.org/archive/

33.

Rupp

A. A.

Templin

(2008). Unique characteristics of diagnostic classification models: A comprehensive review of the current state-of-the-art. Measurement: Interdisciplinary Research & Perspective, 6, 219–262. https://doi.org/10.1080/15366360802490866

34.

Rupp

A. A.

Templin

J. L.

Henson

R. A.

(2010). Diagnostic measurement: Theory, methods and applications. Guilford Press.

35.

Sen

Cohen

A. S.

(2021). Sample size requirements for applying diagnostic classification models. Frontiers in Psychology, 11, 1–16. https://doi.org/10.3389/fpsyg.2020.621251

36.

Sessoms

Henson

R. A.

(2018). Applications of diagnostic classification models: A literature review and critical commentary. Measurement: Interdisciplinary Research and Perspectives, 16, 1–17. https://doi.org/10.1080/15366367.2018.1435104

37.

Sideridis

G. D.

Tsaousis

Al-Harbi

(2021). Identifying student subgroups as a function of school level attributes: A multilevel latent class analysis. Frontiers in Psychology, 12, 1–15. https://doi.org/10.3389/fpsyg.2021.624221

38.

Snow

R. E.

(1991). Aptitude-treatment interaction as a framework for research on individual differences in psychotherapy. Journal of Consulting and Clinical Psychology, 59(2), 205–216. https://doi.org/10.1037/0022-006X.59.2.205

39.

Templin

Hoffman

(2013). Obtaining diagnostic classification model estimates using Mplus. Educational Measurement: Issues and Practice, 32, 37–50. https://doi.org/10.1111/emip.12010

40.

Toker

Green

(2021). A comparison of latent class analysis and the mixture Rasch model using 8th grade mathematics data in the fourth International Mathematics and Science Study (TIMSS-2011). International Journal of Assessment Tools in Education, 8(4), 959–974. https://doi.org/10.21449/ijate.1024251

41.

Tuominen

Niemivirta

Lonka

Salmela-Aro

(2020). Motivation across a transition: Changes in achievement goal orientations and academic well-being from elementary to secondary school. Learning and Individual Differences, 79, 101854. https://doi.org/10.1016/j.lindif.2020.101854

42.

Vermunt

J. K.

Magidson

(2005). Latent GOLD 4.0 user’s guide. Statistical Innovations Inc.

43.

von Davier

(2008). A general diagnostic model applied to language testing data. The British Journal of Mathematical and Statistical Psychology, 61, 287–307. https://doi.org/10.1348/000711007X193957

44.

von Davier

Lee

Y.-S.

(2019). Handbook of diagnostic classification models: Models and model extensions, applications, software packages. Springer Nature. https://doi.org/10.1007/978-3-030-05584-4

45.

Wang

(2021). Using penalized EM algorithm to infer learning trajectories in latent transition CDM. Psychometrika, 86, 167–189. https://doi.org/10.1007/s11336-020-09742-1

46.

White

Murphy

T. B.

(2014). BayesLCA: An R package for Bayesian latent class. Journal of Statistical Software, 61, 1–28. https://doi.org/10.1080/01621459.2016.1231612

47.

(2017). Identifiability of restricted latent class models with binary responses. The Annals of Statistics, 45(2), 675–707. https://doi.org/10.1214/16-AOS1464

48.

Shang

(2018). Identifying latent structures in restricted latent class models. Journal of the American Statistical Association, 113(523), 1284–1295. https://doi.org/10.1080/01621459.2017.1340889

49.

Zhang

(2016). Identifiability of diagnostic classification models. Psychometrika, 81(3), 625–649. https://doi.org/10.1007/s11336-015-9471-z

50.

Yamaguchi

Martinez

A. J.

(in press). Variational Bayes inference for hidden Markov diagnostic classification models. British Journal of Mathematical and Statistical Psychology. https://doi.org/10.1111/bmsp.12308

51.

Yamaguchi

Okada

(2018). Comparison among cognitive diagnostic models for the TIMSS 2007 fourth grade mathematics assessment. PLoS One, 13, e0188691. https://doi.org/10.1371/journal.pone.0188691

52.

Yamaguchi

Okada

(2020a). Variational Bayes inference algorithm for the saturated diagnostic classification model. Psychometrika, 85, 973–995. https://doi.org/10.1007/s11336-020-09739-w

53.

Yamaguchi

Okada

(2020b). Variational Bayes inference for the DINA model. Journal of Educational and Behavioral Statistics, 45(5), 569–597. https://doi.org/10.3102/1076998620911934

54.

Yamaguchi

Templin

(2022a). Direct estimation of diagnostic classification model attribute mastery profiles via a collapsed Gibbs sampling algorithm. Psychometrika, 87(4), 1390–1421.https://doi.org/10.1007/s00357-021-09392-7

55.

Yamaguchi

Templin

(2022b). A Gibbs sampling algorithm with monotonicity constraints for diagnostic classification models. Journal of Classification, 39, 24–54. https://doi.org/10.1007/s00357-021-09392-7

56.

Zhan

Jiao

Liao

(2018). Cognitive diagnosis modelling incorporating item response times. British Journal of Mathematical and Statistical Psychology, 71, 262–286. https://doi.org/10.1111/bmsp.12114

57.

Zhan

Man

Wind

S. A.

Malone

(2022). Cognitive diagnosis modeling incorporating response times and fixation counts: Providing comprehensive feedback and accurate diagnosis. Journal of Educational and Behavioral Statistics, 47(6), 736–776. https://doi.org/10.3102/107699862211110852

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.15 MB