Efficient Metropolis-Hastings Robbins-Monro Algorithm for High-Dimensional Diagnostic Classification Models

Abstract

The expectation-maximization (EM) algorithm is a commonly used technique for the parameter estimation of the diagnostic classification models (DCMs) with a prespecified Q-matrix; however, it requires O(2^K) calculations in its expectation-step, which significantly slows down the computation when the number of attributes, K, is large. This study proposes an efficient Metropolis-Hastings Robbins-Monro (eMHRM) algorithm, needing only O(K + 1) calculations in the Monte Carlo expectation step. Furthermore, the item parameters and structural parameters are approximated via the Robbins-Monro algorithm, which does not require time-consuming nonlinear optimization procedures. A series of simulation studies were conducted to compare the eMHRM with the EM and a Metropolis-Hastings (MH) algorithm regarding the parameter recovery and execution time. The outcomes presented in this article reveal that the eMHRM is much more computationally efficient than the EM and MH, and it tends to produce better estimates than the EM when K is large, suggesting that the eMHRM is a promising parameter estimation method for high-dimensional DCMs.

Keywords

diagnostic classification models stochastic parameter estimation high dimensionality

In recent decades, many various diagnostic classification models (DCMs)—alias cognitive diagnostic models (CDMs)—for item responses have been developed, including, for example, the log-linear cognitive diagnosis model (LCDM; Henson et al., 2009) and its submodels. The successful applications of the DCMs for educational measurement involve extracting information about students’ latent attributes (or called latent knowledge status) to evaluate their skill profiles (e.g., mastery/non-mastery). Beyond the educational environment, the DCMs have been applied in many large-scale tests, such as the Reading and Listening data of the TOEFL^® Internet-based test (IBT), pathological gambling, and the Third International Math and Science Study (TIMSS) (Sessoms & Henson, 2018).

The marginal maximum likelihood with the expectation-maximization (EM) algorithm (Dempster et al., 1977) is a standard parameter estimation method for DCMs, which is accessible from open-source packages such as GDINA (Ma & de la Torre, 2020). The EM algorithm guarantees monotonic climbing to a local minimum of the likelihood function. However, a potential problem with the EM is that it converges very slowly when there is a large amount of “missing data” (i.e., latent attributes). Another approach is to estimate the parameters using the Markov chain Monte Carlo (MCMC) algorithm, such as the Metropolis-Hastings (MH) sampling or the Gibbs sampling. However, the MCMC algorithms could be computationally intensive, especially for many of the attributes (Wang et al., 2021).

Large-scale tests with a large number of latent attributes have been noticed. Tatsuoka et al. (2004) analyzed the 1999 TIMSS mathematics test items from eighth-grade students in 20 countries through the rule-space method (Tatsuoka, 1983). A team of domain experts was invited to identify the probable attributes of 163 items. They developed a Q-matrix of test items that include six content attributes (e.g., geometry and algebra), 10 process attributes (e.g., rule and logic), 11 skill attributes (e.g., unit conversion and proportional reasoning), for a total of 27 attributes. Statistically, the Q-matrix contains element q_ik for item i (i = 1, …, I) and attribute k (k = 1, …K), which can be either 0 or 1, and it is used to indicate whether attribute k is required to solve item i. If attribute k is required, q_ik is 1; otherwise, it is 0. Because the attributes are binary (either mastery or non-mastery), the total number of attribute profiles is 2²⁷. In another example, Lee et al. (2011) analyzed the 2007 TIMSS fourth-grade mathematics test of 25 items and developed a Q-matrix with eight “number” attributes, four “geometric shapes and measures” attributes, and three “data display” attributes, resulting in 2¹⁵ attribute profiles. With such a large number of attributes, it is expected that the EM would be slow.

The present study proposes an efficient Metropolis-Hastings Robbins-Monro (eMHRM) algorithm for saturated DCMs (SDCMs), with a prespecified Q-matrix, to address the problem of many latent attributes. Cai (2010) has shown that the MHRM is more computationally efficient than the EM for multidimensional item response theory models. The numerical integration (for continuous latent traits) or summation (for discrete latent attributes) in the expectation-step of the EM method is not required. Also, it does not need high-intensity sampling for all model parameters used in MCMC methods. The MHRM algorithm uses the MCMC method, partly to sample latent traits, and the maximum likelihood method, to estimate item parameters efficiently. Meanwhile, the MHRM algorithm uses the Robbins-Monro (RM) algorithm (Robbins & Monro, 1951) to gradually average out the Monte Carlo errors caused by Monte Carlo sampling. In this article, the eMHRM advances the conventional MHRM to handle the high-dimensional SDCMs with many attributes.

This article is outlined as follows. The SDCMs are briefly introduced, followed by a review of the EM algorithm. The eMHRM is then developed in detail. Simulation studies are conducted to examine the performance between algorithms in terms of parameter recovery and execution time. Finally, the concluding remarks are given for the eMHRM and its future extensions.

Saturated Diagnostic Classification Models

In recent decades, SDCMs have been actively developed, such as the LCDM. In the present study, the following SDCM is adopted (Yamaguchi & Templin, 2021) as the basis for developing the eMHRM:

\Pr (Y_{n i} = y_{n i} | η_{i}, α_{n}, q_{i}) = \prod_{h = 1}^{2^{\sum_{k = 1}^{K} q_{i k}}} {[η_{i h}^{y_{n i}} {(1 - η_{i h})}^{1 - y_{n i}}]}^{I (d_{i h} = α_{n})}

(1)

This SDCM is formulated where y_ni is the nth respondent’s observed item response to item i (y_ni = 1 for correct response and y_ni = 0 for incorrect response), η_ih denotes the hth probability parameter for the correct response for item i, α_n is the nth respondent’s latent attribute vector indicating his/her mastering status (0: not mastered vs. 1: mastered), q_i represents the ith row vector of a Q-matrix for item i indicating which attributes are required, d_ih is a binary design vector for η_ih, and I(·) denotes the indicator function which equals one if d_ih = α_n and zero otherwise. Taking an item with q = (1, 0, 1) for example (Yamaguchi & Templin, 2021), the $2^{\sum_{k = 1^{q i k}}^{K}} = 2^{2} = 4$ , which means there are four probability parameters for this item: η₁, η₂, η₃, and η₄. The binary design vectors for h = 1, 2, 3, 4 are d₁ = (0, *, 0), d₂ = (0, *, 1), d₃ = (1, *, 0), and d₄ = (1, *, 1), where “*” means the second attribute is irrelevant. The resulting design matrix, called G-matrix (Yamaguchi & Okada, 2020), is shown in Table A1 of Appendix A. The G-matrix denotes a one-to-one mapping between the four probability parameters and a person’s latent attribute status. The outcome in the G-matrix will be one if all the relevant attributes (i.e., first and third attributes) in d equal the corresponding latent attributes; otherwise, it will be zero. The G-matrix is mainly used for convenience of calculation. The major advantage of the SDCM is that the item response probability is characterized as a model parameter. Hence, making use of this feature can easily build up the posterior distributions of model parameters with appropriate conjugate prior distributions. In addition, notice that once the η estimates are obtained, it is feasible to derive the corresponding one-to-one parameter estimates of other SDCMs such as the LCDM (see Yamaguchi & Templin, 2021). Take the LCDM for example, the η₁, η₂, η₃, and η₄ in Table A1 are equated as

\begin{array}{l} η_{1} = {[1 + \exp (- λ_{0})]}^{- 1}, \\ η_{2} = {[1 + \exp (- (λ_{0} + λ_{3}))]}^{- 1}, \\ η_{3} = {[1 + \exp (- (λ_{0} + λ_{1}))]}^{- 1}, \\ η_{4} = {[1 + \exp (- (λ_{0} + λ_{1} + λ_{3} + λ_{13}))]}^{- 1} \end{array}

(2)

where λ₀, λ₁, λ₃, and λ₁₃ are the parameters of the LCDM. Both have the same number of parameters. For the generalized DINA models (de la Torre, 2011), the identity link can be used instead of the logit link adopted above.

The SDCM needs monotonicity constraints on η for each item to identify the model parameters (Yamaguchi & Templin, 2021). For instance, the constraint table for an item with q = (1, 1, 1) is shown in Table B1 of Appendix B. The η₁ must be within 0 and the minimum value among η₂, η₃, and η₅, because the latter three probability parameters require an additional attribute to be mastered. The same logic is applied to the rest of the probability parameters.

The SDCM framework also includes other simpler models as a special case. For instance, the Deterministic Input Noisy “And” gate model (DINA; Macready & Dayton, 1977) can be obtained by re-designing the G-matrix. As shown in Table A2 of Appendix A, there are only two probability parameters for each item: η₁ and η₂. The η₁ corresponds to the guessing parameter, and the 1–η₂ corresponds to the slipping parameter in the DINA model. For illustration purposes, the SDCM and the DINA model will be used as an example in the subsequent simulation studies.

Expectation-Maximization Algorithm

The marginal log-likelihood function of the item and structural parameters is proportional to the following:

l (η, π; y, Q) \propto \sum_{n = 1}^{N} \log [\sum_{c = 1}^{C} π_{c} \prod_{i = 1}^{I} \Pr (y_{n i} | α_{c}, η_{i}, q_{i})]

(3)

Here, N is the number of examinees, I is the number of items, $π_{c}$ is the prior probability of the $α_{c}$ for class c, $\sum_{c = 1}^{c} π_{c} = 1, 0 \leq π_{c} \leq 1$ , and C = 2^K. Instead of the optimizing equation (3), the EM algorithm is often used in practice. The expectation step (E-step) of the EM proposes to optimize a surrogate function, called Q-function, which is proportionally lower bound on the marginal log-likelihood function (Dempster et al., 1977):

Q (ζ | ζ^{(t)}) = \sum_{n = 1}^{N} \sum_{c = 1}^{C} \Pr (α_{c} | y_{n}, ζ^{(t)}, q) [\log π_{c} + \sum_{i = 1}^{I} l o g P r (y_{n i} | α_{c}, η_{i}, q_{i})]

(4)

Here, ζ ϵ (η, π), $\Pr (α_{c} | y_{n}, ζ^{(t)}, q)$ denotes the posterior distribution of α_c given y_n, q, and $ζ^{(t)}$ , which contains the provisional parameter estimates from the previous iteration. Next, the maximization step (M-step) proceeds on maximizing Q(ζ|ζ^$(t)$):

ζ^{(t + 1)} = \underset{ζ}{a r g m a x} Q (ζ | ζ^{(t)})

(5)

The E-step and M-step are carried out alternatively until the sequence ζ^(t) converges to the marginal maximum likelihood estimates. However, the posterior distribution of α_c requires O (2^K) calculations, which becomes intensive in computation when K is large (Wang et al., 2021).

Efficient Metropolis-Hastings Robbins-Monro Algorithm

Yamaguchi and Templin (2021) developed a Gibbs sampler for simulating the attribute pattern for each person. However, that sampler needs to evaluate the multinomial distribution over all attribute patterns (2^K) in the parameter space, which immediately becomes computationally intensive when K is large. In this study, a modified Metropolis-Hastings sampler (Madigan et al., 1995)—for generating the attribute patterns—is employed for reducing the computational burden, which has been used for Bayesian variable selection and DINA models (Liu et al., 2020). The item and structural parameters are then simulated, conditioned on the sampled α values rather than using the mode of their posterior distributions, because the posteriors for item parameters are truncated, whose mode could be outside the parameter space range. The two steps above are iterated until the trajectory of the parameters is stable, which is the very spirit of the data augmentation algorithm (Tanner & Wong, 1987). Finally, the Monte Carlo errors are averaged out through the Robbins-Monro algorithm to obtain stationary point estimates.

The following steps are iterated until specified convergence criteria are satisfied:

\begin{array}{l} Step 1 : α^{(t + 1)} \sim \Pr (α | y, Q, η^{(t)}, π^{(t)}), \\ Step 2 : η^{(t + 1)} \sim \Pr (η | y, Q, α^{(t + 1)}, π^{(t)}) using the RM algorithm, \\ Step 3 : π^{(t + 1)} \sim \Pr (π | y, Q, α^{(t + 1)}, η^{(t + 1)}) using the RM algorithm, \end{array}

(6)

where

t \geq 0

indexes the iterations. The details for each step are described below.

Step 1

Sampling latent attributes via a modified Metropolis–Hastings sampler

The posterior distribution of the current α_n is proportional to the likelihood function multiplying the prior (class) distribution of the current α_n, given by the following:

\Pr (α_{n} = α_{c} | y_{n}, ζ, Q) \propto π_{c}^{I (α_{n} = α_{c})} \prod_{i = 1}^{I} \prod_{h = 1}^{2^{\sum_{k = 1}^{K} q_{i k}}} {[η_{i h}^{y_{n i}} {(1 - η_{i h})}^{1 - y_{n i}}]}^{I (d_{i h} = α_{n})}

(7)

The modified MH sampler suggests the following proposal distribution for a binary attribute:

p [α^{(new)} | f (α^{(t)})] = f {(α^{(t)})}^{α^{(new)}} {[1 - f (α^{(t)})]}^{1 - α^{(new)}}

(8)

which is a Bernoulli distribution with parameter f(α^(t)). The modified MH sampler suggests using f(α^(t)) = 1 – α^(t), which means it always proposes changing the current state at each iteration for each attribute k (i.e., α^(new) = 1 if α^(t) = 0 or α^(new) = 0 if α^(t) = 1). The next step is to determine whether the proposed value of α^(new) is accepted. The MH algorithm constructs an acceptance probability:

\begin{array}{l} a = \min {\frac{\Pr (α^{(new)} | y, ζ, Q) p (α^{(t)} | α^{(new)})}{\Pr (α^{(t)} | y, ζ, Q) p (α^{(new)} | α^{(t)})}, 1} \\ = \min {\frac{\Pr (α^{(new)} | y, ζ, Q)}{\Pr (α^{(t)} | y, ζ, Q)}, 1} \end{array}

(9)

In the above equation, the two proposal distributions are canceled out due to their symmetry, and α^(new) is accepted (i.e., α^(t+1) = α^(new)) when a ≥ U ∼ Uniform(0, 1), otherwise it is rejected (i.e., α^(t+1) = α^(t)).

Notice that the MH only requires O(K+1) operations in coding for k = 1, …, K, because the kth numerator of the acceptance probability (after updated by the acceptance rule) can be re-used for the (k+1)th denominator, except that for the very first iteration the denominator and numerator need to be computed, respectively.

Step 2

Approximating η via the RM algorithm

The posterior distribution for the η_ih is given by

\begin{array}{l} \Pr (η_{i h} | y_{i}, rest) \propto \Pr (y_{i} | η_{i h}, q_{i}, α^{(t + 1)}) \Pr (η_{i h} | a_{i h}, b_{i h}) T (L_{η_{i h}}, U_{η_{i h}}) \\ \propto Beta (a_{i h} + H_{i h}, b_{i h} + T_{i h}) T (L_{η_{i h}}, U_{η_{i h}}), \end{array}

(10)

where

\begin{array}{l} H_{i h} = \sum_{n = 1}^{N} y_{n i} I (d_{i h} = α_{n}), \\ T_{i h} = \sum_{n = 1}^{N} (1 - y_{n i}) I (d_{i h} = α_{n}), \end{array}

(11)

a_ih and, b_ih are the shape parameters of the prior beta distribution for η_ih, and

T (L_{η^{i h}}, U_{η^{i h}})

denotes that the sample space of η_ih is confined between

L_{η^{i h}}

and

U_{η^{i h}}

(see Appendix B). The monotonicity constraints are thus respected. However, the mode (i.e., maximum a posterior (MAP)) of the truncated beta distribution could be outside the lower bound or upper bound during iterations. Instead of finding the mode, the samples of the η are drawn from the truncated beta distribution. This procedure is known as data augmentation (Tanner & Wong, 1987), in that the η estimates are sampled rather than the MAP. Consequently, the computations of gradients, Hessian matrix, and iterative optimization are not required, which can largely reduce the computational burden. In order to restrain boundary estimates (i.e., η_ih = 1 or 0), weak prior information (e.g., a_ih = b_ih = 1.05) is used in this study.

The estimate for the η_ih is error-corrupted due to the noisy input information of the α^(t+1). We utilize the RM algorithm to average out the noisy errors (Cai, 2010; Liu et al., 2020; Robbins & Monro, 1951). Let $g (\hat{η}) = \hat{η} - μ_{\hat{η}}$ and ${\hat{η}}^{(t + 1)} = μ_{\hat{η}} + e^{(t)}$ , where the $g (\hat{η})$ is a function of the difference between $\hat{η}$ (i.e., the noisy estimate) and $μ_{\hat{η}}$ (i.e., $E (\hat{η})$ ), ${\hat{η}}^{(t + 1)} = M_{t}^{- 1} \sum_{m = 1}^{M_{t}} {\hat{η}}^{(t + 1, m)}$ , $M_{t}$ denotes the number of Monte Carlo size at iteration t, and $e^{(t)}$ represents the current measurement error due to the sampling of α^(t+1). Hence, the problem of estimating the average of the $\hat{η}$ (i.e., $μ_{\hat{η}} = E (\hat{η})$ ) is considered as the root-finding of the $g (\hat{η})$ . The estimate is then updated as a weighted mean of the previous estimates and the current one via the RM algorithm, yielding the following:

η_{i h}^{(t + 1)} = η_{i h}^{(t)} + g_{t} (M_{t}^{- 1} \sum_{m = 1}^{M_{t}} {\hat{η}}_{i h}^{(t + 1, m)} - η_{i h}^{(t)}) .

(12)

Here, g_t is the gain sequence defined as

g_{t} \in (0, 1], \sum_{t = 1}^{\infty} g_{t} = \infty, and \sum_{t = 1}^{\infty} g_{t}^{2} < \infty .

(13)

In this case, g_t has the form $b / {(t + 1)}^{w}$ , where $b > 0$ and $1 / 2 < w < 1$ , which yields a g_t, satisfying the conditions in Equation (13).

Step 3

Approximating π via the RM algorithm

The posterior distribution for the π is a Dirichlet distribution, given by the following:

\begin{array}{l} \Pr (π | α^{(t + 1)}, δ, rest) \propto \Pr (α^{(t + 1)} | π) \Pr (π | δ) \\ \propto \prod_{c = 1}^{C} π_{c}^{M_{c} + δ_{c} - 1} \end{array}

(14)

where δ denotes the hyper-parameter of a prior Dirichlet distribution. The MAP estimator for the π might be used. However, the MAP is not stable due to extremely sparse latent attributes when K becomes large while the sample size is relatively small. The samples of the π are generated from the posterior Dirichlet distribution rather than using the MAP. Likewise, the noisy estimate of the

{\hat{π}}_{c}

is updated via the RM algorithm, yielding the following:

π_{c}^{(t + 1)} = π_{c}^{(t)} + g_{t} (M_{t}^{- 1} \sum_{m = 1}^{M_{t}} {\hat{π}}_{c}^{(t + 1, m)} - π_{c}^{(t)})

(15)

Three steps of the eMHRM algorithm are iterated until user-specified termination criteria are satisfied (Yang & Cai, 2014). The first stage sets b = 1 and w = 0 (i.e., g_t = 1), which is considered a “burn-in” period because they are noisy because of the rough starting values. Using w = 0 can make a large step, which usually brings the eMHRM in the vicinity of a solution with a few iterations. The second stage continues with w = 0, where the average of the iterates are used as starting values for the next stage. For the final stage, the iterative averaging procedure with $\frac{1}{2} < w < 1$ was shown as promising for finite-sample problems (Polyak & Juditsky, 1992; Ruppert, 1991), which has also been successfully applied to the parameter estimation of latent variable models (Zhang & Chen, 2020). Specifically, ${\bar{θ}}_{t} = {(t - s)}^{- 1} \sum_{h = s + 1}^{t} {\hat{θ}}_{h}$ shows this, where each of the ${\hat{θ}}_{h}$ is computed as in the typical RM algorithm and s denotes the “burn-in” length. Theoretical results for the iterative averaging show that $\sqrt{t} ({\bar{θ}}_{t} - θ^{*})$ has an asymptotical normal distribution with mean zero and a covariance matrix as small as possible (Polyak & Juditsky, 1992; Ruppert, 1991), where $θ^{*}$ is the solution. The iterative averaging procedure is adopted in the final stage due to its stability and efficiency (Polyak & Juditsky, 1992; Ruppert, 1991; Zhang & Chen, 2020). The final stage continues with w = 1 (Cai, 2010), and the average of the iterates are used as the final estimates.

Regarding standard error and latent attribute estimation, as the eMHRM is a stochastic variant of the EM algorithm, the observed information matrix of the marginal likelihood function can be obtained by Louis (1982) identity (Cai, 2010). For the η and π, their posterior distributions are respectively beta distribution and Dirichlet distribution; therefore, the Hessian matrix and gradient vector are readily available. First, draw a set of Monte Carlo samples for α, given the estimates of η and π. Next, the Louis’s identity is computed based on the α samples (Yang & Cai, 2014). Moreover, the α samples are used to approximate the expected a posteriori (EAP) for the point estimates of α.

Metropolis-Hastings Algorithm

In addition to the EM and the eMHRM, a new MH algorithm is proposed here as a second baseline for comparison purposes. Yamaguchi and Templin (2021) proposed a Gibbs sampler for the parameter estimation of the SDCMs. However, their Gibbs sampler needs to evaluate the posterior distribution of α over 2^K latent classes, requiring O(2^K) calculations. The computation can be burdensome when K is large. A variant of the Gibbs sampler is proposed here. The idea is to replace the Gibbs sampler with the modified MH sampler, which only requires O(K + 1) calculations and is proven more efficient in mutation rates (i.e., the rate of changing α = 0 to α = 1, and vice versa) than the Gibbs sampler (Liu, 1996). For sampling the item parameters and structural parameters, the same posterior distributions described above are used (Yamaguchi & Templin, 2021).

Simulation Study

This simulation study focused on comparing the parameter recovery and execution time for the EM, MH, and eMHRM on the SDCM and DINA model. Once the item and structural parameter estimates are obtained, the latent attributes can be estimated by the EAP estimator. The number of attributes was set to 3, 7, and 15, respectively. The sample sizes were 500, 2,000, and 10,000. The true values of η were generated by $0.2 + (0.8 - 0.2) \times K_{i}^{*} / K_{i}$ (Xu & Shang, 2018), where $K_{i}^{*}$ is the number of attributes out of the total required number of attributes $K_{j}$ for ith item. The entries in α were generated:

α_{k} = {\begin{cases} 1 \\ 0 \end{cases} \begin{array}{l} if θ_{k} \geq 0 \\ otherwise \end{array}, k = 1, . . ., K

(16)

where θ is a continuous variable vector distributed as a multivariate normal distribution with zero means and equal correlations (i.e., ρ = 0, 0.3, and 0.7). The same Q-matrices from Table D1 of Appendix D in Wang et al. (2021) were used, containing 40 items for K = 3, 7, and 15, where the parameters can be identified (Chen et al., 2020; Gu & Xu, 2020; Wang et al., 2021; Xu, 2017; Xu & Zhang, 2016).

Regarding the settings of the eMHRM, the maximum numbers of iterations (200, 200, 900) were for the burn-in stage, averaging stage, and RM stage, respectively. Due to the use of stochastic imputation for latent variables, it was found that the trajectory usually jumps very quickly to the vicinity of the solution (Cai, 2010); therefore, 200 for the burn-in period would be sufficient. The successive 200 iterates were used as the starting values for the final stage. In the final stage, the stopping criterion is set equal to the maximum absolute differences of the averaged item estimates across successive three iterations less than .0001. The M was set equal to 2 for K = 3 and 7, which could stabilize the parameter estimation (Liu et al., 2020), whereas M was set equal to one for K = 15, because we found that the η estimates could be sensitive to M for K = 15. The reason might be that the latent attribute space is highly sparse, and M is often a relatively small value; thus, averaging over the limited number of attribute samples could be sensitive.

For the MH algorithm, the burn-in period was 2,000, and the sample collection period was 3,000, where the length of the burn-in period has been verified sufficient for the SDCM and the DINA model (Wang et al., 2021). The settings of the EM were set by default in the GDINA package. The monotonicity constraints were imposed on the three algorithms.

After algorithm convergence, the person’s attributes were estimated based on the model’s parameter estimates. The EAP estimates were used, where the number of sampling attributes for the eMHRM was 200, which is sufficient, as shown in the subsequent Results section.

The simulation replication number was R = 100 for all conditions, except that for the EM with K = 15 the R was 30 because the EM is too slow and consumed a significant amount of computer memory. Regarding the starting values, the η were randomly generated based on the monotonicity constraints. Also, the starting values of the latent attributes, α, were randomly generated. For the structural parameters, π, the starting values are all set equal to 1/C. The hyper-parameter of the Dirichlet prior, δ_c, for π was set equal to 1, which was used in the MH and eMHRM.

The parameter estimators were assessed by the bias and root-mean-square error (RMSE). Specifically, the bias and the root mean square error (RMSE) of an estimator $\hat{ζ}$ were computed as $R^{- 1} \sum_{r = 1}^{R} (\hat{ζ} - ζ)$ and $\sqrt{R^{- 1} \sum_{r = 1}^{R} {(\hat{ζ} - ζ)}^{2}}$ , respectively, where ζ is the true parameter. For comparison purposes, the mean absolute biases and mean RMSE across η were calculated for each item. For attribute recovery, the two widely used indices, attribute-wise agreement rate (AAR) and pattern-wise agreement rate (PAR), were applied to compute the accuracy of the classification (Wang et al., 2021). Regarding execution time, the absolute execution time was compared among algorithms. The computer’s CPU for the simulation study was 3.7 GHz AMD 3970X. While running the simulation, the maximum CPU speed for one replication was around 4.0 GHz.

Results

For the DINA model, the parameter recoveries of item parameters, η, and structural parameters, π, were shown in Appendices C and D, respectively; the results for the SDCM were shown in Appendices E and F. For the DINA model with K = 3 and 7 (see Appendix C), it is evident that the three algorithms yielded similar mean absolute biases and mean RMSEs for η, regardless of the sample size and correlations. However, for K = 15, the EM could not obtain stable estimates in most cases. In contrast, the eMHRM and MH tended to yield more accurate estimates than the EM. An exception is observed in Figure C9 for K = 15, N = 10,000, and .7 correlation that the EM tended to obtain slightly lower mean absolute biases and RMSEs in most items. The reason might be that the MH and eMHRM used uniform Dirichlet prior, δ = 1, for π, which is not consistent with the true non-uniform distribution of π (i.e., generated based on .7 correlation). Appendix D presents the results of class probability recovery for K = 3, 7, and 15. Overall, the eMHRM and MH yielded similar, or slightly better, results than the EM in some cases (e.g., K = 7 and 15). For .3 and .7 correlations, both ends of the classes (all zero attributes and all one attributes) had slight biases and RMSEs for the three algorithms, which are consistent with the findings in the literature (e.g., Yamaguchi & Okada, 2020). For the SDCM, most of the recovery patterns were similar to those for the DINA model (see Appendices E and F), except that the EM had slightly higher RMSEs for items measuring three attributes for K = 3 (e.g., 12 to 16 items and 28 to 32 items in Figure E1 for N = 500) and K = 7 (e.g., 20 to 24 items and 37 to 40 items in Figure E1 for N = 500). For K = 15, the EM could not perfectly recover the parameters of the SDCM.

Regarding the attribute recovery, Table 1 shows the AAR and PAR for the DINA model and SDCM with zero correlation. For K = 3, the three algorithms obtained similar attribute recovery rates across sample sizes and generating models. For K = 7 and 15, the eMHRM and MH tended to yield better attribute recovery rates than the EM. This might be attributed to the fact that the EM had worse item and structural parameter recovery. For .3 and .7 correlations, the result patterns were similar (see Appendix G), except that the EM obtained slightly better recovery for N = 10,000, K = 15, and .7 correlation, because it had slightly better item and structural parameter estimates, as explained above.

Table 1.

Attribute Recovery for DINA Model and SDCM for Zero Correlation.

(N, K)	DINA			SDCM
AAR	EM	MH	eMHRM	EM	MH	eMHRM
(500, 3)	.975	.975	.975	.974	.974	.974
(500, 7)	.918	.924	.925	.896	.905	.905
(500, 15)	.811	.838	.837	.809	.835	.834
(2000, 3)	.975	.975	.975	.975	.975	.975
(2000, 7)	.926	.927	.927	.908	.909	.909
(2000, 15)	.805	.840	.840	.805	.838	.838
(10000, 3)	.975	.975	.975	.976	.976	.976
(10000, 7)	.927	.927	.927	.911	.911	.911
(10000, 15)	.802	.841	.841	.803	.842	.842
PAR
(500, 3)	.930	.930	.930	.925	.927	.927
(500, 7)	.573	.604	.607	.466	.501	.503
(500, 15)	.049	.074	.074	.044	.066	.066
(2000, 3)	.929	.929	.929	.930	.930	.930
(2000, 7)	.612	.614	.615	.511	.515	.516
(2000, 15)	.039	.077	.076	.038	.070	.070
(10000, 3)	.928	.928	.928	.931	.931	.931
(10000, 7)	.615	.615	.615	.526	.526	.526
(10000, 15)	.037	.077	.077	.038	.075	.075

Note. AAR is for attribute-wise agreement rate and PAR is for pattern-wise agreement rate.

Table 2 shows the average computational time in seconds for the three algorithms when the correlations between attributes were zero. The execution time patterns for .3 and .7 correlations were similar, thus omitted. For K = 3, the EM was slightly faster than the eMHRM, followed by the MH, which is expected because the number of classes was only eight. For K = 7, the eMHRM was only faster than the EM for N = 500 and the SDCM. Since K was 15, it is clear that the eMHRM can achieve a much faster execution time than the EM; for instance, it only took around 6 seconds for N = 500 and around 37 seconds for N = 10,000. It is also worth noting that in the case of K=15 the MH also ran faster than the EM.

Table 2.

Comparison of Average CPU Time (in seconds).

Cor = 0	DINA			SDCM
(N, K)	EM	MH	eMHRM	EM	MH	eMHRM
(500, 3)	0.09	4.22	0.86	0.91	6.08	1.24
(500, 7)	0.57	6.17	1.17	9.20	7.70	1.57
(500, 15)	378.31	61.60	5.39	312.14	70.91	5.71
(2000, 3)	0.28	13.26	2.42	0.95	14.31	2.94
(2000, 7)	1.14	19.35	3.57	4.40	20.64	4.05
(2000, 15)	2228.40	111.24	8.58	1553.70	126.08	8.79
(10000, 3)	1.28	68.49	11.96	1.98	68.01	13.08
(10000, 7)	5.86	104.76	17.11	10.47	106.24	18.14
(10000, 15)	42,791.00	475.33	36.82	15,934.00	496.73	31.96

Note. “Cor” is for correlation between attributes, “N” is for sample size, and “K” is for the number of attributes.

Table 3 shows the comparison of average iteration counts and average time of an iteration (in seconds). Recall that the eMHRM had 200 burn-in iterations and 200 averaging iterations. The eMHRM overall can terminate within 50 iterations after the burn-in stage and averaging stage. In most cases, the EM had lower average iteration counts; however, it cost more time per iteration than the eMHRM in all the simulation conditions, especially for K = 15. The results suggest that the eMHRM was more efficient in computation than the EM and had a better convergence speed.

Table 3.

Comparison of average iteration counts and average time of an iteration (in seconds).

Cor = 0	Average Iteration Count				Average Time of an Iteration
Cor = 0	DINA		SDCM		DINA		SDCM
(N, K)	EM	eMHRM	EM	eMHRM	EM	eMHRM	EM	eMHRM
(500, 3)	6.87	432.91	10.40	445.72	0.014	0.002	0.088	0.003
(500, 7)	34.10	435.22	96.80	448.09	0.017	0.003	0.096	0.003
(500, 15)	136.73	447.43	111.90	463.66	2.806	0.012	2.845	0.012
(2000, 3)	6.45	421.00	8.09	430.24	0.043	0.006	0.118	0.007
(2000, 7)	18.84	422.02	32.71	431.32	0.061	0.008	0.135	0.009
(2000, 15)	214.87	429.40	150.63	439.69	10.518	0.020	10.564	0.020
(10000, 3)	6.06	412.78	7.24	417.55	0.212	0.029	0.274	0.031
(10000, 7)	16.34	413.35	26.26	418.24	0.360	0.041	0.399	0.043
(10000, 15)	770.00	417.13	298.57	423.07	56.101	0.088	54.337	0.076

Note. “Cor” is for correlation between attributes, “N” is for sample size, and “K” is for the number of attributes. For the eMHRM, the maximum numbers of iterations (200, 200, 900) are for the burn-in stage, averaging stage, and Robbins-Monro stage, respectively.

Concluding Remarks

The current study proposes an eMHRM algorithm for high-dimensional DCMs with many attributes. The eMHRM replaces the time-consuming E-step of the EM by a stochastic approximation via the modified MH sampler, which simply needs O(K + 1) calculations rather than O(2^K) in the EM. For the structural parameters, π, the eMHRM requires at least one computation of O(2^K), which is faster than the EM which often needs several iterations (including log-likelihood, gradients, and Hessian matrix) to find the optimal estimates. Considering the sample size and the item size, both need O(I) for the item size, but the eMHRM has a complexity about O(NM) and the EM has O(N2^K) for sample size. Furthermore, the item and structural parameters are approximated via the RM algorithm, which also largely reduces computation time since iterative Newton-type optimization is not required (e.g., no need to compute gradients and Hessian matrix). Also, the inequality constraints on η can be fully respected, thanks to the truncated posterior beta distribution. Thus, time-consuming inequality optimization is not needed. More importantly, the eMHRM is customized for the SDCM, which does not need to specify the logistic or probit link function. In contrast, past studies using the Gibbs sampling are restricted to the probit link for constructing exact conditional distributions (e.g., Wang et al., 2021).

In addition to the eMHRM, a new MH algorithm for high-dimensional DCMs is proposed for comparison purposes, a variant of the Gibbs algorithm (Yamaguchi & Templin, 2021). The significant advantage of the MH algorithm is that it only requires O(K + 1) computation in sampling attributes. In contrast, the Gibbs algorithm needs to evaluate 2^K multinomial distribution, which can be slow in practice. The MH algorithm performed well in all simulation conditions regarding the parameter recovery; although not as fast as the eMHRM, its speed might be acceptable in practice for K = 15 (i.e., around 8 minutes).

Regarding computational efficiency, the results presented in Table 2 show that the eMHRM was faster than the MH for a large K (e.g., K = 15), not to mention the EM. Considering the results of attribute recovery (Table 1) and time elapse (Tables 2 and 3), the EM can be adopted for the DINA model with K = 3 as it is faster and yields identical attribute recovery. For K = 7 and K = 15, the EM yields worse attribute recovery than the eMHRM and the MH, though faster in some cases for K = 7; thus it is suggested to adopt the eMHRM for good attribute recovery and fast execution for K = 7 and K = 15.

Two algorithms not included in the present study, the sequential Gibbs sampler (Wang et al., 2021) and the variational EM (Yamaguchi & Okada, 2020), were recently proposed to tackle the parameter estimation of high-dimensional SDCMs with a prespecified Q-matrix. For the SDCM with K = 15, N = 2,000, and I = 40, the eMHRM took 8.79 seconds (under maximum CPU speed of 4.0 GHz), much more efficient than the sequential Gibbs sampler, which requires 106.66 seconds (under maximum CPU speed of 4.7 GHz) (see Table 9 in Wang et al., 2021). Comparing the eMHRM with the variational EM, the variational EM was reported as taking 1000 seconds in N = 5000 and K = 10 conditions, whereas the eMHRM took 31.96 seconds in N = 10,000 and K = 15 conditions. The eMHRM could be faster than the variational EM because the latter must evaluate 2^K classes of the posteriori of the attributes (Yamaguchi & Templin, 2021), which is time-consuming for K = 15.

The parameter recovery performance of the EM was not stable for K = 7, and in most cases, it became worse for K = 15. A possible explanation for the poor performance might be that the EM was trapped in the local maximum of the log-likelihood function, which has been reported by Jiang and Ma (2018). One solution might be using global optimization methods, such as the differential evolution algorithm, to alleviate the local-maximum problem (Jiang & Ma, 2018). However, the execution time would take much longer than the report in Table 2, which is not desirable in practice. This issue is beyond the scope of the current study, so it is left for further research.

The current eMHRM could be extended to more complex issues. For example, estimating the item and structural parameters together with the Q-matrix has been an important issue when some elements of the Q-matrix are not correctly specified by domain experts. As the eMHRM is very efficient in computation, extending the eMHRM for the Q-matrix estimation of the high-dimensional SDCMs is under study for future research. In addition, using the number of Monte Carlo size, M = 2 for K = 3 and 7, and M = 1 for K = 15 seem sufficient for our simulation studies considered; however, more investigation for different combinations of K and M is required to examine the sensitivity of parameter estimates.

Supplemental Material

Supplemental Material - Efficient Metropolis-Hastings Robbins-Monro Algorithm for High-Dimensional Diagnostic Classification Models

Supplemental Material for Efficient Metropolis-Hastings Robbins-Monro Algorithm for High-Dimensional Diagnostic Classification Models by Chen-Wei Liu in Applied Psychological Measurement

Footnotes

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, authorship, and/or publication of this article: The first author acknowledges the grant support from the Taiwan Ministry of Science and Technology, Grant Number MOST 110-2410-H-003-054-MY2.

ORCID iD

Chen-Wei Liu

Supplemental Material

Supplemental material for this article is available online.

References

Cai

(2010). High-dimensional exploratory item factor analysis by a Metropolis–Hastings Robbins–Monro algorithm. Psychometrika, 75(1), 33–57. https://doi.org/10.1007/s11336-009-9136-x

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

de la Torre

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

Dempster

A. P.

Laird

N. M.

Rubin

D. B.

(1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 39(1), 1–38. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x

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

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(2), 191–210. https://doi.org/10.1007/s11336-008-9089-5

Jiang

(2018). Integrating differential evolution optimization to cognitive diagnostic model estimation. Frontiers in Psychology, 9, 2142. https://doi.org/10.3389/fpsyg.2018.02142

Lee

Y. S.

Park

Y. S.

Taylan

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

Liu

C. W.

Andersson

Skrondal

(2020). A constrained Metropolis–Hastings Robbins–Monro algorithm for Q matrix estimation in DINA models. Psychometrika, 85(2), 322–357. https://doi.org/10.1007/s11336-020-09707-4

10.

Liu

J. S.

(1996). Peskun’s theorem and a modified discrete-state Gibbs sampler. Biometrika, 83(3), 681–682. https://doi.org/10.1093/biomet/83.3.681

11.

Louis

T. A.

(1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 44(2), 226–233. https://doi.org/10.1111/j.2517-6161.1982.tb01203.x

12.

de la Torre

(2020). GDINA: An R package for cognitive diagnosis modeling. Journal of Statistical Software, 93(14), 1–26. https://doi.org/10.18637/jss.v093.i14

13.

Macready

G. B.

Dayton

C. M.

(1977). The use of probabilistic models in the assessment of mastery. Journal of Educational Statistics, 2(2), 99–120. https://doi.org/10.3102/10769986002002099

14.

Madigan

York

Allard

(1995). Bayesian graphical models for discrete data. International Statistical Review/Revue Internationale de Statistique, 63(2), 215–232. https://doi.org/10.2307/1403615

15.

Polyak

B. T.

Juditsky

A. B.

(1992). Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4), 838–855. https://doi.org/10.1137/0330046

16.

Robbins

Monro

(1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22(3), 400–407. https://doi.org/10.1214/aoms/1177729586

17.

Ruppert

(1991). Stochastic approximation. In Ghosh

B. K.

Sen

P. K.

(Eds.), Handbook in sequential analysis (pp. 503–529). Marcel Dekker.

18.

Sessoms

Henson

R. A.

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

19.

Tanner

M. A.

Wong

W. H.

(1987). The calculation of posterior distributions by data augmentation. Journal of the American statistical Association, 82(398), 528–540. https://doi.org/10.2307/2289460

20.

Tatsuoka

K. K.

(1983). Rule space: An approach for dealing with misconceptions based on item response theory. Journal of Educational Measurement, 20(4), 345–354. https://doi.org/10.1111/j.1745-3984.1983.tb00212.x

21.

Tatsuoka

K. K.

Corter

J. E.

Tatsuoka

(2004). Patterns of diagnosed mathematical content and process skills in TIMSS-R across a sample of 20 countries. American Educational Research Journal, 41(4), 901–926. https://doi.org/10.3102/00028312041004901

22.

Wang

Shi

Zhang

(2021). Sequential Gibbs sampling algorithm for cognitive diagnosis models with many attributes. Multivariate Behavioral Research, 23, 1–37. https://doi.org/10.1080/00273171.2021.1896352

23.

(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

24.

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

25.

Zhang

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

26.

Yamaguchi

Okada

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

27.

Yamaguchi

Templin

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

28.

Yang

J. S.

Cai

(2014). Estimation of contextual effects through nonlinear multilevel latent variable modeling with a Metropolis–Hastings Robbins–Monro algorithm. Journal of Educational and Behavioral Statistics, 39(6), 550–582. https://doi.org/10.76998614559972

29.

Zhang

Chen

(2020). Computation for latent variable model estimation: A unified stochastic proximal framework. Psychometrika, Online First. https://doi.org/10.1007/s11336-022-09863-9

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

6.84 MB