Modeling Skipped Items in Cognitive Diagnostic Assessments

Abstract

In cognitive diagnosis assessments, examinees frequently skip some items due to various reasons, such as time constraints, lack of confidence, or the perception of item difficulty. These skipped items are usually due to lacking specific cognitive attributes or knowledge structures. However, many research studies conventionally rely on complete-case analysis, potentially compromising the validity of results by ignoring skipped items. This omission can lead to invalid inferences about the attributes profiles of examinees and introduce bias into item parameters. This study aims to develop statistical models for effectively handling missing data resulting from skipped items. Specifically, we employ an item response theory model for missing indicators and utilize the deterministic inputs, noisy “and” gate (DINA) model to describe cognitive item responses. Furthermore, we introduce a higher-order structure to describe the correlation between the higher-order ability parameters and skipping propensity parameters and the correlation of the item parameters. The proposed new model fills the gaps in missing data handling, providing a more precise evaluation for the examinees who skip items and yielding more accurate research results. This enhances our understanding of cognitive processes and strengthens support for educational policies, aligning them better with student needs and optimizing decision-making.

Plain language summary

We propose a missing data model for skipped items in cognitive diagnostic assessments. The significance of our findings is illustrated with an application to a real data set.

Introduction

In cognitive diagnostic assessments, researchers typically utilize responses from multiple items to gain a comprehensive understanding of individuals’ cognitive attributes or knowledge structures. However, a common challenge emerges when examinees choose to respond to some items while skipping others, resulting in missing data at item level. This behavior can be attributed to various reasons, including perceived difficulty of specific items, allocating excessive time to one item leading to the omission of others, or simply skipping some items. Ignoring such missing data can lead to information loss and potentially bias the interpretations. It is important to note that many research investigations commonly employ an analytical approach that relies solely on available data, evaluating examinees based exclusively on their observed response data. This approach may potentially provide incomplete diagnostic feedback to individuals with missingness in their data, and importantly, may lead to biased conclusions when comparing them to examinees with complete datasets.

The process of making a cognitive diagnosis has received extensive attention in psychology and educational assessment recently. A great deal of cognitive diagnosis models (CDMs) has been developed, for instance, the deterministic inputs, noisy “and” gate (DINA) model (Junker & Sijtsma, 2001), the log-linear cognitive diagnostic (LCDM) model (Henson et al., 2009), and the general DINA model (de la Torre, 2009; de la Torre & Douglas, 2004; Junker & Sijtsma, 2001). De la Torre and Douglas (2004) proposed a higher-order DINA (HO-DINA) model to connect latent attributes through higher-order ability.

In practice, cognitive diagnosis assessments could investigate examinees’ skills and knowledge, which provide fine-grained diagnosis and feedback information, such as attribute mastery information and cognitive processes. For example, even if the scores of two examinees are the same, their knowledge structures and learning paths may be different. Thus, teachers’ instructional approaches should also be distinct. Consequently, examinees with the same score tend to have distinct cognitive structures and cognitive processes. In this case, if examinees’ strengths and weaknesses can be identified, much more meaningful and much richer information about the individuals’ attribute profiles can be provided. Meanwhile, such information can also be used to adapt instruction to the examinees and determine their skills and knowledge that need to be developed further.

Cognitive diagnosis assessments (de la Torre, 2011) are effective tools for evaluating examinees’ mastery of skills or attributes. However, in practical test scenarios, some examinees often leave a few items blank for various reasons. For example, they may perceive some items as too difficult, they may have insufficient time to answer some questions, or their test-taking motivation may decline over the course of the assessments. Usually, the skipped/omitted items occur within a block of responses. According to the Programme for the International Assessment of Adult Competencies (PIAAC) in 2012, skipped item rates varied significantly; in Korea, only 2% of items in the numeracy domain were skipped, while in Chile, 25.9% of items in the literacy domain were left unanswered (OECD, 2013). This highlights the widespread issue of missing data in educational assessment and psychological measurement (Rose et al., 2017; Yuan et al., 2018). Consequently, if these missing responses are not handled properly, they can pose a significant threat to the validity of the assessments (Pohl et al., 2014; Rose et al., 2015).

To address the missing data issue more effectively, Rubin (1976), Little and Rubin (2002) proposed three mechanisms: not missing at random (NMAR), missing at random (MAR), and missing completely at random (MCAR). Let Y be a response matrix with the elements of observed response data ${\mathbf{Y}}_{\mathrm{observe}}$ and missing data ${\mathbf{Y}}_{\mathrm{missing}}$ , and let $\mathbf{D}$ be the corresponding missing indicator matrix. Consequently, MCAR is described as $\mathrm{P}(\mathbf{D}|{\mathbf{Y}}_{\mathrm{observe}},{\mathbf{Y}}_{\mathrm{missing}})=\mathrm{P}(\mathbf{D})$ , MAR is expressed as $\mathrm{P}(\mathbf{D}|{\mathbf{Y}}_{\mathrm{observe}},{\mathbf{Y}}_{\mathrm{missing}})=\mathrm{P}(\mathbf{D}|{\mathbf{Y}}_{\mathrm{observe}})$ , and NMAR denotes that the missing data are correlated with both ${\mathbf{Y}}_{\mathrm{observe}}$ and ${\mathbf{Y}}_{\mathrm{missing}}$ .

In the framework of IRT, the issue of missing data has garnered substantial attention, leading to the development of various methodologies to address item nonresponse based on different missing mechanisms (Dai, 2021; Du et al., 2022; Pohl & Becker, 2020; Sperrin & Martin, 2020; Xiao & Bulut, 2020). The contributions include the multidimensional IRT (MIRT) models explored by scholars (Glas & Pimentel, 2008; Grund et al., 2021;Holman & Glas, 2005; Moustaki & Knott, 2000; O’Muircheartaigh & Moustaki, 1999), as well as latent regression models (LRM) proposed by Rose et al. (2010). The literature has also witnessed sophisticated integrations, such as the combination of MIRT with LRM as presented by Rose et al. (2017), and innovative models like the IRTree proposed by Debeer et al. (2017). A recent development, as emphasized by Pohl et al. (2019) and Ulitzsch et al. (2020a, 2020b), focuses on incorporating response time to model item nonresponses (Guo et al., 2022; Lu et al., 2023; Lu & Wang, 2020).

Moving into the realm of cognitive diagnosis assessments (Ma, 2022; Templin & Henson, 2006), which intricately analyze examinees’ skills and knowledge, the need to address missing data becomes increasingly urgent. While several methods have emerged within the cognitive diagnosis framework, most remain relatively straightforward. Ömür Sünbül (2018) investigated four methods to address missing data in the DINA model under the MCAR and MAR mechanisms. These methods include treating missing data as incorrect, person mean imputation, two-way imputation, and expectation-maximization imputation. While these methods are conceptually clear, their effectiveness tends to diminish as the proportion of missing items increases. Specifically, when dealing with non-ignorable missing data, these methods can introduce biases and errors into statistical inference. Following this, Dai et al. (2018) examined the influence of missing responses on two Q matrix validation methods. However, they did not delve into the impact of missingness on specific CDMs. Shan and Wang (2020) proposed a joint CDM model that addresses item-level missing data by introducing the categorical missingness propensity of examinees, but they did not elaborate on this potential categorization. Dai and Svetina Valdivia (2022) assessed the performance of different strategies for managing missing data, including the full information maximum likelihood estimation in CDMs. However, their focus remained on methods to treat missing responses in the CDMs, without constructing a dedicated model for these missing responses. More recently, Liang et al. (2022) introduced a joint cognitive model for not-reached missingness. Yet, when it comes to skipped items, to our knowledge, no model has been developed within the cognitive diagnoses framework to particularly address them.

Recognizing the absence of such models, this paper introduces a novel missing data model for skipped items within the CDM framework. Our aim extends beyond merely addressing this gap; we aim to present a model that offers a comprehensive understanding, departing from the more simplistic approaches employed previously. In essence, our work not only assesses skipped items but also delivers a nuanced analysis, representing a significant advancement. To capture the underlying processes of the cognitive diagnosis, we bring together research on modeling cognitive item responses and research on missing data models in a joint modeling framework. This framework is characterized by a bivariate normal distribution between examinees’ higher-order ability parameters and skipping propensity parameters, and a multivariate normal distribution of the item intercept parameters, interaction parameters, and item skipping difficult parameters. The parameters of the model were estimated by the Bayesian Markov chain Monte Carlo (MCMC) method. The results of simulations showed that the proposed model improved attribute profile correct classification rates and produced precise estimation of item parameters when the missing mechanism is nonignorable.

The remainder of this paper is structured as follows. First, a missing indicator model is presented to model the skipped item responses. Second, the higher-order DINA model is used to model the response data. The missing model and response model are linked by the correlation of higher-level ability parameters and the skipping propensity parameters, and the correlation of item parameters. Third, the MCMC algorithm was used to estimate the model parameters and two model selection criteria were given to evaluate the model fit. Fourth, simulation studies are conducted to evaluate the performance of the proposed model under different conditions. Finally, concluding remarks and some suggestions for future research are provided.

Model Construct

Let $\mathit{Y}$ be a two-dimensional data matrix, if we observe the response data for the $i$ th examinee to the $j$ th item, $Y_{\mathrm{ij}}$ corresponds to the observation $y_{\mathrm{ij}}$ . Otherwise, it is missing data, and we use the sign “s” to express the missing data and the relevant parameters.

Modeling Skipped Items

The missing data mechanism of skipped items is modeled using an additional manifest or latent variable, which represents the examinees’ tendency to skip items (Holman & Glas, 2005; Rose et al., 2017). The skipping indicator $D_{\mathrm{ij}}$ contains information on whether examinee $i$ generated a response to item $j$ , with 1 indicating a skipped response and 0 representing an observed response. The missing indicator matrix D is defined as:

$$D_{\mathrm{ij}}=\{\begin{array}{c}\multicolumn{1}{c}{0,\mathrm{if}{y}_{\mathrm{ij}}\mathrm{was}\mathrm{observed},}\\ \multicolumn{1}{c}{1,\mathrm{if}{y}_{\mathrm{ij}}\mathrm{was}\mathrm{not}\mathrm{observed}.}\end{array}$$

(1)

The Rasch model (Rasch, 1966) is employed to express the probability of skipping items. That is, the probability of item skipping is modeled as a function of an examinee’s skipping propensity ${\theta }_i^s$ and the $j$ th item skipping difficulty ${\beta }_j^s$ . The probability of item skipping is modeled as:

$$p(s_{\mathrm{ij}}=1|{\theta }_i^s,{\beta }_j^s)=\frac{\exp ({\theta }_i^s-{\beta }_j^s)}{1+\exp ({\theta }_i^s-{\beta }_j^s)}$$

(2)

HO-DINA Model

The DINA model describes the probability of item response as a function of latent attributes, and the probability of the $i$ th examinee responding to item $j$ correctly is expressed as follows:

$$p(Y_{\mathrm{ij}}=1)=g_j+(1-s_j-g_j)\underset{k=1}{\overset{K}{\Pi }}{{\alpha }_{\mathrm{ik}}}^{q_{\mathrm{jk}}},$$

(3)

where $g_j$ and $s_j$ are the guessing and slipping probabilities of the $j$ th item, respectively; $1-s_j-g_j=\mathrm{ID}{I}_j$ is the $j$ th item discrimination index (De La Torre, 2008); and ${\alpha }_{\mathrm{ik}}$ describes the attribute mastery statue for examinee $i$ ,that is, the value of ${\alpha }_{\mathrm{ik}}$ is 1 when examinee $i$ masters attribute $k$ , otherwise, the value of ${\alpha }_{\mathrm{ik}}$ is 0. The Q matrix (Tatsuoka, 1983) is a $J\times K$ matrix with element $q_{\mathrm{jk}}$ ; when the attribute $k$ is required for answering the $j$ th item correctly, the value of $q_{\mathrm{jk}}$ is 1, otherwise it is 0.

Equation 3 can be re-parameterized by the logit scale as the re-parameterized DINA model (DeCarlo, 2011):

$${\beta }_j=\mathrm{logit}(g_j),$$

(4)

$${\delta }_j=\mathrm{logit}(1-s_j)-\mathrm{logit}\left(g_j\right).$$

(5)

In addition, $\mathrm{logit}(x)=\log (\frac{x}{1-x})$ , thus, Equation 3 can be expressed as:

$$\mathrm{logit}(P(y_{\mathrm{ij}}=1))={\beta }_j+{\delta }_j\underset{k=1}{\overset{K}{\Pi }}{\alpha }_{\mathrm{ik}}^{q_{\mathrm{jk}}},$$

(6)

where ${\beta }_j$ and ${\delta }_j$ are the item intercept and interaction parameters, respectively. The higher-order structure is flexible: on one hand, it can produce more accurate higher-order ability and attribute structure estimates; in addition, it also reduces the number of model parameters. On the other hand, the attributes are usually correlated, and it would be more beneficial to formulate a higher-order structure to link attributes (de la Torre & Douglas (2004); Zhan et al., 2018), which is expressed as:

$$\mathrm{logit}(P({\alpha }_{\mathrm{ik}}=1))={{\theta }_i}^h{\gamma }_k-{\lambda }_k,$$

(7)

where $P({\alpha }_{\mathrm{ik}}=1)$ represents the probability that the $i$ th examinee possesses the $k$ th attribute, ${\theta }_i^h$ is the higher-order ability for examinee $i,$ and ${\gamma }_k$ and ${\lambda }_k$ are the slope and intercept parameters for attribute $k$ , respectively. The slope parameter ${\gamma }_k$ is positive, indicating that higher values of an ability ${\theta }_i^h$ correspond to better knowledge mastery of the attribute by the examinee.

The Hierarchical Model

Following the hierarchical modeling framework, item parameters ${\beta }_j$ and ${\delta }_j$ of the DINA model and item skipping difficult parameter ${\beta }_j^s$ are assumed to follow a multivariate normal distribution with the mean vector and covariance matrix, which is expressed as follows:

$$\left(\begin{array}{c}{\beta }_j\hfill \\ \multicolumn{1}{c}{{\delta }_j}\\ \multicolumn{1}{c}{{\beta }_j^s}\end{array}\right)~N\left((\begin{array}{c}{\mu }_{\beta }\hfill \\ \multicolumn{1}{c}{{\mu }_{\delta }}\\ \multicolumn{1}{c}{{\mu }_{{\beta }^s}}\end{array},{\mathbf{\Sigma }}_{\mathbf{I}})\right),{\mathbf{\Sigma }}_{\mathbf{I}}=\left(\begin{array}{ccc}{\sigma }_{\beta }^2\hfill & {\sigma }_{\beta \delta }\hfill & {\sigma }_{\beta {\beta }^s}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\beta \delta }}& {\sigma }_{\delta }^2\hfill & {\sigma }_{\delta {\beta }^s}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\beta {\beta }^s}}& {\sigma }_{\delta {\beta }^s}\hfill & {\sigma }_{{\beta }^s}^2\hfill \end{array}\right)$$

(8)

In addition, the joint distribution of the person parameters is assumed to follow a bivariate normal distribution with mean vector ${\mathbf{\mu }}_{\mathrm{P}}$ and covariance matrix ${\mathbf{\Sigma }}_{\mathrm{P}}$ :

$${\mathbf{\mu }}_{\mathrm{P}}=({\mu }_{{\theta }^h},{\mu }_{{\theta }^s})\prime ,{\mathbf{\Sigma }}_{\mathrm{P}}=\left(\begin{array}{cc}{\sigma }_{{\theta }^h}^2\hfill & {\sigma }_{{\theta }^h{\theta }^s}\hfill \\ \multicolumn{1}{c}{{\sigma }_{{\theta }^h{\theta }^s}}& {\sigma }_{{\theta }^s}^2\hfill \end{array}\right).$$

(9)

Model Identifications

To solve the identifications of the proposed model, we have made the following restrictions. On the one hand, to eliminate the trade-off between ability ${\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}$ and skipping threshold parameter ${\mathrm{\beta }}_{\mathrm{j}}^{\mathrm{s}}$ , we set ${\mathrm{\mu }}_{{\mathrm{\theta }}^{\mathrm{s}}}=0$ . Similarly, to eliminate the trade-off between the higher-order ability person parameter ${\mathrm{\theta }}^{\mathrm{h}}$ and the attribute intercept ${\mathrm{\lambda }}_{\mathrm{k}}$ , we set ${\mathrm{\mu }}_{{\mathrm{\theta }}^{\mathrm{h}}}=0$ . On the other hand, ${\mathrm{\sigma }}_{{\mathrm{\theta }}^{\mathrm{h}}}$ =1 was fixed to eliminate the scale trade-off between ${\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}$ and ${\mathrm{\gamma }}_{\mathrm{k}}$ . These restrictions were based on the method to fix the mean of ability to 0 and the variance of ability to 1 (Fox, 2010; Lord & Novick, 1968).

Model Selection

We used two common Bayesian model evaluation criteria: the deviance information criteria (DIC; Spiegelhalter et al., 2002) and the logarithm of the pseudo marginal likelihood (LPML, Geisser & Eddy, 1979; Ibrahim et al., 2001); They can be employed for comparing the missing mechanism models through MCMC sampling. The $\mathrm{DIC}$ is given by:

$$\mathrm{De}{\mathrm{v}}_{(\mathit{Y},\mathit{D}|\mathbf{\Omega })}=-2\log L(\mathit{Y},\mathit{D}|\mathbf{\Omega })$$

$$=-2\sum _{i=1}^N\sum _{j=1}^J\sum _{k=1}^K[(Y_{\mathrm{ij}}=s)\log (P(Y_{\mathrm{ij}}=s))$$

$$+(Y_{\mathrm{ij}}=1)\log ((1-P(Y_{\mathrm{ij}}=s))P(Y_{\mathrm{ij}}=1))$$

$$+(Y_{\mathrm{ij}}=0)\log ((1-P(Y_{\mathrm{ij}}=s))P(Y_{\mathrm{ij}}=0))].$$

(10)

Based on the posterior distribution of $\mathbf{Dev}(\mathit{Y},\mathit{D}|\mathit{\Omega })$ , the $\mathrm{DIC}$ is defined as:

$$\mathbf{DIC}=\overline{\mathbf{Dev}}+{\mathit{p}}_{\mathit{D}}=\overline{\mathbf{Dev}}+(\overline{\mathbf{Dev}}-\widehat{\mathbf{Dev}}),$$

(11)

where $\overline{\mathbf{Dev}}=\mathit{E}(\mathbf{Dev}(\mathit{Y},\mathit{D},\mathit{\Omega })|\mathit{Y},\mathit{D})\cong \frac{1}{\mathit{R}}\sum _{\mathit{r}=1}^{\mathit{R}}\mathbf{Dev}(\mathit{Y},\mathit{D},{\mathit{\Omega }}^{\mathit{r}})$ , which is the posterior mean deviance; a Bayesian measure of fit, $\mathit{r}=1,...,\mathit{R}$ denotes the $r$ th iteration of the algorithm; and $\widehat{\mathbf{Dev}}=\mathbf{Dev}(\mathit{Y},\mathit{D},\overline{\Omega })$ , which is the effective number of parameters, is a Bayesian measure of complexity, with $\overline{\Omega }=\mathit{E}(\mathit{\Omega }|\mathit{Y},\mathit{D})\cong \frac{1}{\mathit{R}}\sum _{\mathit{r}=1}^{\mathit{R}}{\mathit{\Omega }}^{\mathit{r}}$ . A smaller $\mathbf{DIC}$ indicates a better model fit.

Another way to compare the two model fits is LPML. The conditional predictive ordinate index (Chen et al., 2000) of the two models is computed. Let ${\mathit{Q}}_{\mathit{ij},\mathit{max}}=\mathit{ma}{\mathit{x}}_{1\le \mathit{r}\le \mathit{R}}\{-\mathit{logf}({\mathit{Y}}_{\mathit{ij}},{\mathit{D}}_{\mathit{ij}}|{\mathit{\Omega }}^{\mathit{r}})\}$ . Then,

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{log}\left(\widehat{\mathit{CP}{\mathit{O}}_{\mathit{ij}}}\right)=-{\mathit{Q}}_{\mathit{ij},\mathit{max}}}\\ \multicolumn{1}{c}{-\mathit{log}\left[\frac{1}{\mathit{R}}\sum _{\mathit{r}=1}^{\mathit{r}}\mathit{exp}\{-\mathit{logf}({\mathit{Y}}_{\mathit{ij}},{\mathit{D}}_{\mathit{ij}}|{\mathit{\Omega }}^{\mathit{r}})-{\mathit{Q}}_{\mathit{ij},\mathit{max}}\}\right].}\end{array}$$

(12)

Note that the maximum value adjustment used in $\mathbf{log}\left(\widehat{\mathit{CP}{\mathit{O}}_{\mathit{ij}}}\right)$ plays an important role in numerical stabilization in computing $\mathit{exp}\{-\mathit{logf}({\mathit{Y}}_{\mathit{ij}},{\mathit{D}}_{\mathit{ij}}|{\mathit{\Omega }}^{\mathit{r}})-{\mathit{Q}}_{\mathit{ij},\mathit{max}}\}$ in Equation 12. The summary statistic for log $\left(\widehat{\mathit{CP}{\mathit{O}}_{\mathit{ij}}}\right)$ is the sum of their logarithms, and the LPML is:

$$\mathbf{LPML}=\sum _{\mathit{i}=1}^{\mathit{N}}\sum _{\mathit{j}=1}^{\mathit{J}}\mathit{log}({\widehat{\mathit{CPO}}}_{\mathit{ij}})$$

(13)

where the model with a larger LPML signifies a better fit to the data.

Simulation Studies

We conducted four simulation studies to evaluate different aspects of the proposed model. Simulation study I was conducted to assess whether the Bayesian MCMC algorithm could successfully recover parameters for the proposed model under different numbers of examinees with fixed numbers of items and attributes. Simulation study II was conducted to investigate the parameter recovery of different numbers of attributes for the same numbers of examinees and items. Simulation study III was intended to show the differences in model parameter estimates between the NMAR and MAR models for different skipping proportions and correlations between person parameters. Simulation study IV aimed to show the performance of the model selection using DIC and LPML for the NMAR and MAR models with different skipping proportions and correlations between person parameters.

We employed the Bayesian MCMC method to estimate the model parameters. The convergence of the Bayesian algorithm was checked by monitoring the trace plots of the parameters for 10,000 iterations. The first 5,000 iterations were set as the burn-in period. There were R = 30 replications for each simulation condition. In terms of convergence, we found that the potential scale reduction factor (PSRF; Brooks & Gelman, 1998) was less than 1.1 for all parameters, which means that the chains converge well. The bias and root mean square error (RMSE) were computed to evaluate the parameter recoveries. To evaluate item parameter recoveries, the bias was averaged across all items. Regarding the classification of individual attributes and profiles, the attribute correct classification rate (ACCR) and the pattern correct classification rate (PCCR) were computed.

Data Generation

In the four simulation studies, the item parameters were sampled from the following distributions:

$$\begin{array}{c}\multicolumn{1}{c}{\left(\begin{array}{c}{\beta }_j\hfill \\ \multicolumn{1}{c}{{\delta }_j}\\ \multicolumn{1}{c}{{\beta }_j^s}\end{array}\right)~\mathrm{MVN}(\left(\begin{array}{c}{\mu }_{\beta }\hfill \\ \multicolumn{1}{c}{{\mu }_{\delta }}\\ \multicolumn{1}{c}{{\mu }_{{\beta }^s}}\end{array}\right),{\mathbf{\Sigma }}_{\mathit{I}}),{\mu }_{\beta }=-2.197,{\mu }_{\delta }=4.394,}\\ \multicolumn{1}{c}{{\mathbf{\Sigma }}_{\mathbf{I}}=\left(\begin{array}{ccc}1\hfill & -0.8\hfill & 0.25\hfill \\ \multicolumn{1}{c}{-0.8}& 1\hfill & -0.15\hfill \\ \multicolumn{1}{c}{0.25}& -0.15\hfill & 0.25\hfill \end{array}\right),}\end{array}$$

which were consistent with Liang et al. (2022). The low and high skipping proportions were varied by setting the value of ${\mu }_{{\beta }^s}$ . That is, the skipping proportion was 12% (low) when ${\mu }_{{\beta }^s}=2.1$ ; it was 25% (high) when ${\mu }_{{\beta }^s}=1.3.$

The intercept parameters of the attributes were set as $\mathit{\lambda }=(-1,-0.5,0,0.5,1),$ and the slope parameters were ${\gamma }_k=1.5$ for all attributes, consistent with Liang et al. (2022). We considered three Q matrices (Shan & Wang, 2020; Xu & Shang, 2018), as expressed in Figure 1.

Figure 1.

$Q_1$ presents a matrix with 30 items and five attributes, $Q_2$ represents a matrix with 20 items and five attributes, and $Q_3$ represents a matrix with 20 items and three attributes.

The person parameters ${\theta }_i^h$ and ${\theta }_i^s$ were simulated from the bivariate normal distribution $\left(\begin{array}{c}{\theta }^h\hfill \\ \multicolumn{1}{c}{{\theta }^s}\end{array}\right)~\mathrm{MVN}(\left(\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right),\left(\begin{array}{cc}1\hfill & {\sigma }_{{\theta }^h{\theta }^s}\hfill \\ \multicolumn{1}{c}{{\sigma }_{{\theta }^h{\theta }^s}}& {\sigma }_{{\theta }^s}^2\hfill \end{array}\right)),\mathrm{where}{\sigma }_{{\theta }^s}^2$ = 0.25. Three levels of the correlation between ${\theta }_i^h$ and ${\theta }_i^s$ (i.e., ${\mathit{\rho }}_{{\theta }_i^h{\theta }_i^s}$ ) were considered:0 (uncorrelated), −.5 (medium), and −.8 (high). Two levels of skipping proportions were 12% (low) and 25% (high).

$$Q_1=\left(\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 1\hfill & 1\hfill \end{array}\right)Q_2=\left(\begin{array}{ccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right)Q_3=\left(\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 1\hfill & 1\hfill \end{array}\right)$$

The Prior and Posterior Distributions

Model parameters were estimated using the Bayesian MCMC algorithm. The joint posterior distributions for the model parameters are obtained through observed data likelihood and prior distributions of model parameters.

Let $\mathit{\Omega }=\{{\theta }_i^h,{\theta }_i^s,{\beta }_j^s,{\alpha }_{\mathrm{ik}},{\beta }_j,{\delta }_j,{\gamma }_k,{\lambda }_k,{\mu }_{\beta },{\mu }_{\delta },{\mu }_{{\beta }^s},{\mathbf{\Sigma }}_{\mathbf{I}},{\sigma }_{{\theta }^h{\theta }^s},{\sigma }_{{\theta }^s}^2\}$ be the parameters of the proposed model. The priors and hyperpriors of the model parameters are presented in Table 1.

Table 1.

The Priors and Hyperpriors of the Model Parameters and Their Values in Simulation Studies.

Parameters	Priors/hyperpriors	Values
${\mu }_{\beta }$	$N({\mu }_{\beta },{\sigma }_{\beta }^2)$	${\mu }_{\beta }=-2.197,{\sigma }_{\beta }^2=2$
${\mu }_{\delta }$	$N({\mu }_{\delta },{\sigma }_{\delta }^2)I({\mu }_{\delta }>0)$	${\mu }_{\delta }=4.394,{\sigma }_{\delta }^2=2$
${\mu }_{\beta }^s$	$N({\mu }_{{\beta }^s},{\sigma }_{{\beta }^s}^2)$	${\mu }_{{\beta }^s}=2.1,{\sigma }_{{\beta }^s}^2=2$
${\lambda }_K$	$N({\mu }_{\lambda },{\sigma }_{\lambda }^2)$	${\mu }_{\lambda }=0,{\sigma }_{\lambda }^2=4$
${\gamma }_K$	$N({\mu }_{\gamma },{\sigma }_{\gamma }^2)I({\gamma }_K>0)$	${\mu }_{\gamma }=0,{\sigma }_{\gamma }^2=4$
${\mathbf{\Sigma }}_I$	$\mathrm{Inv}-\mathrm{Wishart}({\mathbf{\Sigma }}_{I0}^{-1},v_{I0})$	${\mathbf{\Sigma }}_{I0}=\left(\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right),v_{I0}=2,k_{I0}=1$
${\sigma }_{{\theta }^h{\theta }^s}$	${\sigma }_{{\theta }^h{\theta }^s}~U(-a,a)$	$a=1$
${\sigma }_{{\theta }^s}^2$	Inv $-\mathrm{\gamma }(v_0,c_0)$	$v_0=2,c_0=2$

The posterior distribution of the proposed model can be obtained as follows:

$$\begin{array}{c}\multicolumn{1}{c}{p(\mathit{\Omega }|\mathit{Y},\mathit{D})=\underset{i=1}{\overset{N}{\Pi }}\underset{j=1}{\overset{J}{\Pi }}\underset{k=1}{\overset{K}{\Pi }}P(Y_{\mathrm{ij}}=s{)}^{Y_{\mathrm{ij}}=s}}\\ \multicolumn{1}{c}{\times {\left\{[1-P(Y_{\mathrm{ij}}=s)]P(Y_{\mathrm{ij}}=1)\right\}}^{Y_{\mathrm{ij}}=1}}\\ \multicolumn{1}{c}{\times {\left\{[1-P(Y_{\mathrm{ij}}=s)]P(Y_{\mathrm{ij}}=0)\right\}}^{Y_{\mathrm{ij}}=0}\times p({\theta }^h,{\theta }^s;{\mathit{\mu }}_{\mathbf{P}},{\mathbf{\Sigma }}_{\mathbf{P}})}\end{array}$$

$$p({\beta }_j,{\delta }_j,{\beta }_j^s;{\mu }_{\beta },{\mu }_{\delta },{\mu }_{{\beta }^s},{\mathbf{\Sigma }}_{\mathbf{I}})p\left({\mathbf{\Sigma }}_{\mathbf{P}}\right)p\left({\mu }_{\beta }\right)p\left({\mu }_{\delta }\right)p({\mu }_{{\beta }^s})p({\mathbf{\Sigma }}_{\mathbf{I}})$$

$$p({\alpha }_{\mathrm{ik}};{\theta }_i^h,{\gamma }_k,{\lambda }_k)p\left({\lambda }_k\right)p\left({\gamma }_k\right).$$

(14)

The details of the MCMC algorithm are given in the Appendix.

Simulation Study I

In simulation study I, the aim was to investigate the impact of the sample size on the performance of the parameters and attribute profiles. The conditions were set as follows: 1,000 examinees responded to 30 items with five attributes, and this was compared to a scenario where 500 examinees answered 30 items with five attributes. Moreover, the correlation between ${\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}$ and ${\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}$ was −.3, and the skipping proportion was low (12%).

Table 2 presents recoveries of the ability parameters and item parameters as well as the attribute parameter estimates. Given the test length and the number of attributes, when the number of examinees increased from 500 to 1,000, the bias and RMSE of the item parameters $\mathit{\beta },\mathit{\delta }$ , and ${\mathit{\beta }}^{\mathit{s}}$ decreased. Specifically, 30 items and five attributes were fixed, the number of examinees increased from 500 to 1,000, and the average RMSE of all the item intercept parameters ranged from 0.233 to 0.191. For the interaction parameters, the average bias of all item interaction parameters ranged from 0.014 to 0.001, and the average RMSE ranged from 0.366 to 0.291. For the skipping item difficulty parameters, the average bias ranged from 0.021 to 0.018, and the average RMSE ranged from 0.133 to 0.099. Moreover, the RMSE of the item mean parameters and item covariance matrix elements decreased noticeably. In addition, the first two columns of Table 4 show the ACCRs and PCCRs for simulation study I. With the test length and number of attributes held constant, an increase in the number of examinees from 500 to 1,000 led to higher ACCRs and PCCRs, owing to the larger sample size. Notably, the PCCRs exhibited a distinct variation, ranging from 0.848 to 0.852. Figure 2 shows the PSRF values of several items with 500 examinees, 30 items, and five attributes. We found that the item intercept parameters $\mathit{\beta }$ , the item interaction parameters $\mathit{\delta },$ and the skipping difficulty parameters ${\mathit{\beta }}^{\mathit{s}}$ converged at 5,000 iterations.

Table 2.

The Parameter Estimates in Simulation Study I.

	N = 1,000		N = 500
	J = 30		J = 30
	K = 5		K = 5
Parameter	Bias	RMSE	Bias	RMSE
$\beta$	0.011	0.191	0.007	0.233
$\delta$	0.001	0.291	0.014	0.366
${\beta }^s$	0.018	0.099	0.021	0.133
${\mu }_{\beta }$	0.036	0.179	0.032	0.186
${\mu }_{\delta }$	−0.129	0.228	−0.118	0.236
${\mu }_{{\beta }^s}$	0.016	0.091	0.021	0.102
${\lambda }_1$	0.012	0.148	−0.014	0.213
${\lambda }_2$	−0.024	0.135	−0.079	0.232
${\lambda }_3$	−0.043	0.128	−0.003	0.195
${\lambda }_4$	−0.023	0.156	−0.077	0.223
${\lambda }_5$	−0.014	0.177	−0.019	0.248
${\gamma }_1$	0.041	0.298	0.036	0.381
${\gamma }_2$	0.014	0.265	0.037	0.378
${\gamma }_3$	−0.001	0.234	0.166	0.516
${\gamma }_4$	0.037	0.297	−0.046	0.416
${\gamma }_5$	0.043	0.279	−0.014	0.399
${\sigma }_{\beta }^2$	−0.129	0.298	−0.148	0.319
${\sigma }_{\delta }^2$	−0.043	0.325	−0.069	0.352
${\sigma }_{{\beta }^s}^2$	−0.043	0.073	−0.033	0.073
${\sigma }_{\beta \delta }$	0.165	0.298	0.207	0.337
${\sigma }_{\beta {\beta }^s}$	−0.077	0.120	−0.077	0.123
${\sigma }_{\delta {\beta }^s}$	0.063	0.112	0.074	0.124
${\sigma }_{{\theta }^h{\theta }^s}$	−0.005	0.044	0.001	0.061
${\sigma }_{{\theta }^s}^2$	0.022	0.042	0.027	0.053
${\theta }^h$	−0.005	0.628	−0.014	0.635
${\theta }^s$	0.002	0.376	0.003	0.372

The bold elements in the table highlight significant differences in the results.

Figure 2.

The trace plots of PSRF values in simulation study I.

In simulation study I, our findings indicate: First, a larger sample size provides a more robust dataset that captures individual responses comprehensively, thereby enhancing the precision of item parameter estimates. For instance, the bias and RMSE for item parameters were significantly reduced when the sample size was increased to 1,000. Second, when focusing on attributes, ACCR reflects the accuracy of classifying each attribute, while PCCR measures the overall accuracy of pattern classification. Our analysis revealed that as the sample size increased, both ACCR and PCCR values improved. For detailed results, please refer to Table 4. These findings emphasize the significance of sample size in improving the precision of item parameter estimation and enhancing the accuracy of attribute classification.

Simulation Study II

The aim of simulation study II was to investigate the parameter recoveries of different attributes under two conditions: (a) 500 examinees, 20 items, five attributes; and (b) 500 examinees, 20 items, three attributes. In both cases, the correlation between ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ and ${\mathit{\theta }}_{\mathit{i}}^{\mathit{s}}$ was −.3, and the skipping proportion was low (12%).

Table 3 shows the bias and RMSE values for the parameters. The ability parameter ${\mathbf{\theta }}^{\mathbf{h}}$ improved when the number of attributes ranged from K = 3 to K = 5. Furthermore, the ACCRs and PCCRs are displayed in the last two columns of Table 4. Under the condition of 500 examinees and 20 items, a transition from five attributes to three attributes resulted in changes in ACCRs and PCCRs. ACCRs were higher with three attributes, averaging at 0.975, compared to an average of 0.950 with five attributes. Moreover, the PCCRs significantly increased when the number of attributes decreased. Specifically, the PCCR under five attributes was 0.791 compared with the PCCR under three attributes was 0.931. Consequently, both ACCRs and PCCRs were notably higher under three attributes compared to five attributes when considering the fixed condition of 500 examinees and 20 items.

Table 3.

The Parameter Estimates in Simulation Study II.

	N = 500		N = 500
	J = 20		J = 20
	K = 5		K = 3
Parameter	Bias	RMSE	Bias	RMSE
$\beta$	0.011	0.291	−0.007	0.259
$\delta$	0.015	0.428	0.046	0.354
${\beta }^s$	0.022	0.145	0.022	0.134
${\mu }_{\beta }$	−0.142	0.275	−0.160	0.276
${\mu }_{\delta }$	0.109	0.289	0.138	0.291
${\mu }_{{\beta }^s}$	−0.111	0.191	−0.112	0.191
${\lambda }_1$	0.047	0.215	−0.025	0.216
${\lambda }_2$	−0.021	0.193	−0.025	0.198
${\lambda }_3$	−0.015	0.198	−0.022	0.202
${\lambda }_4$	−0.011	0.202	—	—
${\lambda }_5$	−0.067	0.257	—	—
${\gamma }_1$	0.032	0.473	0.074	0.664
${\gamma }_2$	0.034	0.398	0.211	0.782
${\gamma }_3$	0.209	0.531	0.116	0.551
${\gamma }_4$	0.035	0.398	—	—
${\gamma }_5$	−0.023	0.386	—	—
${\sigma }_{\beta }^2$	−0.028	0.424	−0.065	0.389
${\sigma }_{\delta }^2$	0.170	0.549	0.162	0.501
${\sigma }_{{\beta }^s}^2$	0.180	0.238	0.177	0.235
${\sigma }_{\beta \delta }$	0.018	0.410	0.038	0.373
${\sigma }_{\beta {\beta }^s}$	0.032	0.194	0.027	0.178
${\sigma }_{\delta {\beta }^s}$	−0.008	0.199	−0.006	0.191
${\sigma }_{{\theta }^h{\theta }^s}$	0.001	0.065	−0.006	0.069
${\sigma }_{{\theta }^s}^2$	0.046	0.071	0.045	0.068
${\theta }^h$	−0.008	0.640	−0.021	0.705
${\theta }^s$	0.010	0.396	0.010	0.395

The bold elements in the table highlight significant differences in the results.

Table 4.

ACCRs and PCCRs in Simulation Studies I and II.

N = 1,000	N = 500	N = 500	N = 500
J = 30	J = 30	J = 20	J = 20
K = 5	K = 5	K = 5	K = 3
0.961	0.960	0.926	0.964
0.956	0.956	0.946	0.982
0.964	0.964	0.968	0.980
0.961	0.962	0.948	—
0.980	0.979	0.964	—
0.852	0.848	0.791	0.931

The bold elements in the table highlight significant differences in the results.

In summary, a larger sample size and longer item length resulted in improved ACCR and PCCR recoveries. Conversely, when the number of examinees or the item length was reduced, both ACCRs and PCCRs decreased. Note that these changes became more apparent when the test length was decreased.

Simulation Study III

Simulation study III aims to investigate the parameter recoveries under different skipping proportions and correlation coefficients. We considered three models: the NMAR, MAR, and HO-DINA models. Two skipping proportions were considered: 12% (low) and 25% (high). In addition, three correlations of ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{h}}$ and ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{s}}$ were 0 (uncorrelated), −0.5 (medium), and −0.8 (high). As a result, there were 2 × 3 simulation conditions. The data were generated from the proposed model under the NMAR mechanism. There were 30 items answered by 500 examinees, and each item associated with five attributes.

For the conditions with a low skipping proportion, Tables 5 and 6 show the bias and RMSE values for the parameters, including person, item, and attribute parameters of the three models, based on the different correlations between ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{s}}$ and ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{h}}$ . When the correlation between ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{s}}$ and ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{h}}$ was 0, the parameter estimates for the three models were quite similar. However, the NMAR model exhibited an advantage as the correlation between ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ and ${\mathit{\theta }}_{\mathit{i}}^{\mathit{s}}$ increased. For instance, the bias and RMSEs of parameters $\mathit{\beta },\mathit{\delta },{\mathit{\mu }}_{\mathit{\beta }},{\mathit{\mu }}_{\mathit{\delta }},\lambda ,\gamma ,{\mathit{\theta }}^{\mathit{h}},\mathbf{and}{\mathit{\theta }}^{\mathit{s}}$ for the NMAR model were consistently the smallest in the NMAR model compared to the MAR and HO-DINA models. As the correlation between ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{h}}$ and ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{s}}$ increased, the performance of the attribute parameters became particularly evident. For example, the RMSE of the attribute parameters $\mathbf{\lambda }\mathbf{and}\mathbf{\gamma }$ in the NMAR model were smaller compared with those of the MAR and HO-DINA models. As expected, the NMAR model demonstrated higher accuracy in terms of parameter estimation. Furthermore, because ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{h}}$ and ${\mathbf{\theta }}_{\mathbf{i}}^{\mathbf{s}}$ in both the MAR and HO-DINA models were uncorrelated, we observed that the parameter estimates in these models were relatively consistent in simulations that did not account for skipped items.

Table 5.

The Bias of Parameter Estimates in the Models With Low Skipping Proportions and Different ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ in Simulation Study III.

Parameter	$\mathrm{\rho }=$ 0			$\mathrm{\rho }=-$ 0.5			$\mathrm{\rho }=-$ 0.8
	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA
$\beta$	0.014	0.011	0.010	0.002	0.167	0.167	0.017	0.178	0.177
$\delta$	0.001	0.005	0.025	0.001	−0.063	−0.046	0.001	−0.061	−0.041
${\beta }^s$	0.019	0.017	—	0.019	0.018	—	0.015	0.013	—
${\mu }_{\beta }$	0.041	0.037	0.037	0.027	0.191	0.192	0.045	0.201	0.203
${\mu }_{\delta }$	−0.133	−0.127	−0.107	−0.132	−0.195	−0.178	−0.134	−0.191	−0.173
${\mu }_{{\beta }^s}$	0.018	0.017	—	0.018	0.017	—	0.015	0.011	—
${\sigma }_{\beta }^2$	−0.113	−0.106	−0.100	−0.131	−0.212	−0.213	−0.142	−0.223	−0.218
${\sigma }_{\delta }^2$	−0.056	−0.052	0.145	−0.155	−0.201	−0.019	−0.097	−0.150	0.029
${\sigma }_{{\beta }^s}^2$	−0.031	−0.031	—	−0.036	−0.036	—	−0.030	−0.030	—
${\sigma }_{\beta \delta }$	0.170	0.164	0.144	0.216	0.294	0.288	0.199	0.287	0.279
${\sigma }_{\beta {\beta }^s}$	−0.075	−0.075	—	−0.079	−0.078	—	−0.078	−0.075	—
${\sigma }_{\delta {\beta }^s}$	0.062	0.062	—	0.078	0.085	—	0.067	0.071	—
${\lambda }_1$	0.042	0.042	0.035	0.051	0.282	0.280	−0.002	0.193	0.196
${\lambda }_2$	−0.031	−0.036	−0.028	−0.022	0.218	0.221	−0.058	0.187	0.197
${\lambda }_3$	0.055	0.055	0.051	−0.021	0.189	0.187	−0.033	0.163	0.164
${\lambda }_4$	0.021	0.017	0.016	−0.013	0.268	0.268	0.012	0.274	0.280
${\lambda }_5$	0.001	0.004	0.008	0.023	0.178	0.174	−0.041	0.092	0.093
${\gamma }_1$	0.070	0.066	0.067	0.092	0.707	0.696	0.154	0.698	0.690
${\gamma }_2$	0.012	0.014	0.027	−0.004	0.488	0.486	0.012	0.562	0.563
${\gamma }_3$	0.110	0.110	0.098	0.119	0.510	0.501	0.075	0.429	0.425
${\gamma }_4$	0.040	0.036	0.028	−0.102	0.344	0.349	−0.012	0.400	0.400
${\gamma }_5$	−0.029	−0.024	−0.012	0.062	0.255	0.248	−0.093	0.070	0.068
${\sigma }_{{\theta }^h{\theta }^s}$	−0.006	—	—	−0.021	—	—	−0.015	—	—
${\sigma }_{{\theta }^s}^2$	0.035	0.031	—	0.043	0.040	—	0.024	0.014	—
${\theta }^h$	0.004	0.003	0.003	−0.008	−0.011	−0.012	−0.008	−0.011	−0.010
${\theta }^s$	0.001	0.001	—	−0.003	−0.001	—	0.001	0.001	—

Note. NMAR = not missing at random model; MAR = missing at random model; HO-DINA = higher-order DINA model.The bold elements in the table highlight significant differences in the results.

Table 6.

The RMSE of Parameter Estimates in the Models With Low Skipping Proportions and Different ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ in Simulation Study III.

Parameter	$\mathrm{\rho }=$ 0			$\mathrm{\rho }=-$ 0.5			$\mathrm{\rho }=-$ 0.8
	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA
$\beta$	0.230	0.232	0.234	0.252	0.385	0.386	0.240	0.372	0.375
$\delta$	0.384	0.385	0.391	0.398	0.461	0.454	0.394	0.446	0.448
${\beta }^s$	0.132	0.132	—	0.144	0.144	—	0.131	0.132	—
${\mu }_{\beta }$	0.188	0.187	0.188	0.182	0.257	0.257	0.184	0.263	0.265
${\mu }_{\delta }$	0.243	0.240	0.245	0.231	0.268	0.267	0.235	0.267	0.269
${\mu }_{{\beta }^s}$	0.099	0.100	—	0.101	0.101	—	0.101	0.102	—
${\sigma }_{\beta }^2$	0.313	0.312	0.316	0.325	0.328	0.328	0.315	0.330	0.331
${\sigma }_{\delta }^2$	0.360	0.358	0.422	0.359	0.364	0.343	0.360	0.360	0.364
${\sigma }_{{\beta }^s}^2$	0.072	0.072	—	0.073	0.073	—	0.072	0.072	—
${\sigma }_{\beta \delta }$	0.324	0.321	0.328	0.341	0.371	0.373	0.333	0.370	0.373
${\sigma }_{\beta {\beta }^s}$	0.124	0.123	—	0.125	0.121	—	0.124	0.120	—
${\sigma }_{\delta {\beta }^s}$	0.118	0.118	—	0.123	0.125	—	0.120	0.121	—
${\lambda }_1$	0.241	0.237	0.239	0.226	0.361	0.358	0.228	0.325	0.325
${\lambda }_2$	0.198	0.199	0.199	0.183	0.293	0.294	0.193	0.293	0.294
${\lambda }_3$	0.191	0.190	0.188	0.202	0.284	0.283	0.187	0.266	0.263
${\lambda }_4$	0.206	0.205	0.206	0.196	0.343	0.343	0.225	0.355	0.358
${\lambda }_5$	0.235	0.236	0.239	0.232	0.285	0.285	0.229	0.247	0.251
${\gamma }_1$	0.429	0.434	0.438	0.396	0.839	0.824	0.445	0.847	0.836
${\gamma }_2$	0.421	0.426	0.434	0.386	0.638	0.639	0.368	0.722	0.716
${\gamma }_3$	0.418	0.417	0.408	0.404	0.668	0.648	0.397	0.603	0.590
${\gamma }_4$	0.404	0.401	0.394	0.412	0.539	0.546	0.391	0.578	0.580
${\gamma }_5$	0.343	0.347	0.351	0.394	0.444	0.441	0.365	0.364	0.369
${\sigma }_{{\theta }^h{\theta }^s}$	0.069	—	—	0.063	—	—	0.054	—	—
${\sigma }_{{\theta }^s}^2$	0.057	0.055	—	0.062	0.061	—	0.051	0.045	—
${\theta }^h$	0.630	0.629	0.629	0.621	0.638	0.638	0.589	0.634	0.633
${\theta }^s$	0.379	0.379	—	0.364	0.381	—	0.329	0.380	—

Note. NMAR = not missing at random model; MAR = missing at random model; HO-DINA = higher-order DINA model.The bold elements in the table highlight significant differences in the results.

For the conditions with a high skipping proportion, as presented in Tables 7 and 8, we observed similar results to those observed in conditions with low skipping proportions. Across both scenarios, the NMAR model consistently outperformed the other two models, particularly in the accuracy of person and item parameter estimates. Attribute parameter estimation also displayed notable differences among the three mechanism models, highlighting the NMAR model’s superior performance.

Table 7.

The Bias of Parameter Estimates in the Models With High Skipping Proportions and Different ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ in Simulation Study III.

Parameter	$\mathrm{\rho }=$ 0			$\mathrm{\rho }=-$ 0.5			$\mathrm{\rho }=-$ 0.8
	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA
$\beta$	0.021	0.020	0.019	0.018	0.282	0.290	0.023	0.291	0.296
$\delta$	0.007	0.008	0.028	−0.017	−0.140	−0.129	−0.019	−0.148	−0.135
${\beta }^s$	0.001	0.0001	—	0.004	0.003	—	0.004	0.004	—
${\mu }_{\beta }$	0.048	0.044	0.047	0.043	0.303	0.313	0.051	0.313	0.318
${\mu }_{\delta }$	−0.126	−0.123	−0.106	−0.150	−0.269	−0.260	−0.155	−0.278	−0.262
${\mu }_{{\beta }^s}$	0.001	−0.002	—	0.003	0.001	—	0.003	0.002	—
${\sigma }_{\beta }^2$	−0.156	−0.152	−0.148	−0.147	−0.208	−0.213	−0.153	−0.217	−0.222
${\sigma }_{\delta }^2$	−0.084	−0.081	0.113	−0.092	−0.116	0.066	−0.089	−0.109	0.071
${\sigma }_{{\beta }^s}^2$	−0.084	−0.084	—	−0.084	−0.084	—	−0.083	−0.084	—
${\sigma }_{\beta \delta }$	0.210	0.207	0.192	0.203	0.272	0.270	0.205	0.272	0.270
${\sigma }_{\beta {\beta }^s}$	−0.074	−0.073	—	−0.071	−0.076	—	−0.074	−0.083	—
${\sigma }_{\delta {\beta }^s}$	0.025	0.025	—	0.017	0.025	—	0.016	0.029	—
${\lambda }_1$	0.038	0.040	0.035	0.017	0.487	0.510	0.023	0.566	0.602
${\lambda }_2$	−0.030	−0.029	−0.024	−0.110	0.451	0.478	−0.055	0.455	0.467
${\lambda }_3$	0.041	0.041	0.037	0.004	0.361	0.362	0.031	0.397	0.389
${\lambda }_4$	0.009	0.006	0.017	0.038	0.485	0.499	0.016	0.440	0.458
${\lambda }_5$	0.011	0.009	0.003	0.004	0.234	0.230	0.003	0.196	0.185
${\gamma }_1$	0.119	0.119	0.104	0.084	1.191	1.230	0.179	1.628	1.717
${\gamma }_2$	−0.016	−0.002	0.0003	0.070	1.285	1.367	0.016	1.159	1.158
${\gamma }_3$	0.101	0.096	0.086	0.080	0.668	0.672	0.0001	0.633	0.627
${\gamma }_4$	0.086	0.081	0.106	0.046	0.733	0.753	0.039	0.689	0.723
${\gamma }_5$	0.014	0.012	0.002	0.010	0.301	0.283	0.018	0.220	0.202
${\sigma }_{{\theta }^h{\theta }^s}$	−0.006	—	—	−0.002	—	—	−0.006	—	—
${\sigma }_{{\theta }^s}^2$	0.024	0.022	—	0.018	0.016	—	0.023	0.022	—
${\theta }^h$	−0.006	−0.005	−0.006	−0.014	−0.023	−0.022	0.002	−0.006	−0.006
${\theta }^s$	−0.002	−0.003	—	−0.001	−0.002	—	0.003	0.003	—

Table 8.

The RMSE of Parameter Estimates in the Models With High Skipping Proportions and Different ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ in Simulation Study III.

Parameter	$\mathrm{\rho }=$ 0			$\mathrm{\rho }=-$ 0.5			$\mathrm{\rho }=-$ 0.8
	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA	NMAR	MAR	HO-DINA
$\beta$	0.272	0.276	0.282	0.247	0.511	0.530	0.267	0.544	0.559
$\delta$	0.428	0.431	0.438	0.390	0.522	0.534	0.394	0.539	0.552
${\beta }^s$	0.108	0.107	—	0.103	0.103	—	0.102	0.102	—
${\mu }_{\beta }$	0.188	0.187	0.190	0.187	0.348	0.357	0.190	0.356	0.362
${\mu }_{\delta }$	0.247	0.246	0.254	0.249	0.332	0.334	0.253	0.339	0.339
${\mu }_{{\beta }^s}$	0.086	0.086	—	0.084	0.084	—	0.085	0.083	—
${\sigma }_{\beta }^2$	0.331	0.335	0.337	0.328	0.326	0.328	0.334	0.337	0.339
${\sigma }_{\delta }^2$	0.395	0.395	0.432	0.357	0.358	0.383	0.377	0.367	0.395
${\sigma }_{{\beta }^s}^2$	0.097	0.097	—	0.097	0.097	—	0.097	0.097	—
${\sigma }_{\beta \delta }$	0.358	0.359	0.364	0.340	0.362	0.368	0.349	0.368	0.375
${\sigma }_{\beta {\beta }^s}$	0.114	0.114	—	0.113	0.113	—	0.114	0.117	—
${\sigma }_{\delta {\beta }^s}$	0.097	0.097	—	0.093	0.093	—	0.093	0.094	—
${\lambda }_1$	0.237	0.234	0.236	0.218	0.551	0.576	0.233	0.641	0.675
${\lambda }_2$	0.196	0.198	0.199	0.218	0.527	0.552	0.195	0.528	0.541
${\lambda }_3$	0.195	0.193	0.193	0.193	0.415	0.419	0.220	0.456	0.448
${\lambda }_4$	0.241	0.238	0.247	0.245	0.551	0.566	0.229	0.505	0.525
${\lambda }_5$	0.260	0.252	0.254	0.239	0.332	0.330	0.256	0.310	0.301
${\gamma }_1$	0.460	0.448	0.445	0.436	1.346	1.404	0.454	1.862	1.922
${\gamma }_2$	0.431	0.445	0.442	0.396	1.457	1.559	0.408	1.376	1.333
${\gamma }_3$	0.462	0.444	0.453	0.431	0.787	0.792	0.354	0.754	0.749
${\gamma }_4$	0.507	0.496	0.521	0.458	0.864	0.882	0.421	0.818	0.856
${\gamma }_5$	0.470	0.462	0.464	0.381	0.448	0.436	0.391	0.430	0.412
${\sigma }_{{\theta }^h{\theta }^s}$	0.052	—	—	0.055	—	—	0.048	—	—
${\sigma }_{{\theta }^s}^2$	0.042	0.041	—	0.040	0.039	—	0.042	0.043	—
${\theta }^h$	0.646	0.646	0.645	0.623	0.654	0.654	0.579	0.650	0.651
${\theta }^s$	0.336	0.336	—	0.321	0.332	—	0.296	0.333	—

The bold elements in the table highlight significant differences in the results.

Figures 3 and 4 show the bias of attribute parameters $\mathrm{\lambda }$ and $\mathrm{\gamma }$ in the NMAR and MAR models across different skipping proportions and correlations between ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ and ${\mathit{\theta }}_{\mathit{i}}^{\mathit{s}}$ . The results indicate that the parameter estimates between the NMAR and MAR models were similar when the correlation between ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ and ${\mathit{\theta }}_{\mathit{i}}^{\mathit{s}}$ was 0. Regardless of correlation level or skipping proportion, the bias of attribute parameters $\mathrm{\lambda }$ and $\mathrm{\gamma }$ remained close to 0 in the NMAR model, indicating accurate and effective estimation of these attributes. In contrast, when the skipping proportion was high, the bias of the attribute parameters $\mathrm{\lambda }$ and $\mathrm{\gamma }$ of the MAR model were significantly larger than those of the NMAR model. In addition, the RMSE of the attribute parameters $\mathbf{\lambda }$ and $\mathbf{\gamma }$ are shown in Figures 5 and 6. Comparing the NMAR and MAR models under these conditions, we observe that the RMSEs of attribute parameters in the NMAR model exhibited significant improvement. The estimates of attribute parameters became notably more precise with increasing correlation level or higher skipping proportions. In contrast, the RMSEs of attribute parameters in the MAR model were exceptionally large.

Figure 3.

Bias of parameter estimates about $\mathrm{\lambda }$ with different skipping proportions and ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ .

Figure 4.

Bias of parameter estimates about $\mathrm{\gamma }$ with different skipping proportions and ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ . .

Figure 5.

RMSE of parameter estimates about $\mathrm{\lambda }$ with different skipping proportions and ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ . .

Figure 6.

RMSE of parameter estimates about $\mathrm{\gamma }$ with different skipping proportions and ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ . .

Figure 7 presents the ACCRs and PCCRs in the NMAR and MAR models under six simulation conditions: low and high skipping proportions, as well as low, medium, high correlations between ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ and ${\mathit{\theta }}_{\mathit{i}}^{\mathit{s}}$ . Detailed results of the ACCRs and PCCRs are provided in Table S1 in the Supplemental Material. The ACCRs and PCCRs in the NMAR model consistently exhibit a noticeable advantage, especially in scenarios with high skipping proportions and high correlations. This indicates that the NMAR model is more effective in recovering attribute patterns when the missing data mechanism is non-ignorable. Conversely, the MAR model struggles to effectively recover attribute patterns in such scenarios.

Figure 7.

The ACCRs and PCCRs for MAR and NMAR models in simulation study III.

In brief, our results consistently indicate that the NMAR model outperforms both the MAR model and the HO-DINA model, regardless of whether the skipping proportion is high or low, or whether the person correlations are strong or weak. Furthermore, across all scenarios, the NMAR model consistently achieves higher ACCRs and PCCRs compared to the other models. This highlights that the NMAR model excels in delivering more precise parameter estimation and attribute profiles.

Simulation Study IV

In simulation study IV, our aim was to conduct the model selection using Bayesian model assessment criteria, namely, DIC and LPML, based on the six conditions from simulation III. We used 500 examinees who answered 30 items with five attributes. The model selection results are shown in Table 9. There were no significant differences in DIC and LPML among the models when ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ was 0. However, when  ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ deviated from 0, the DIC of the NMAR model was consistently smaller, and the LPML was relatively higher. This suggests that the DIC and LPML criteria effectively facilitated the accurate selection of the true model.

Table 9.

Model Selections Under Different ${\mathrm{\rho }}_{{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{h}}{\mathrm{\theta }}_{\mathrm{i}}^{\mathrm{s}}}$ and Skipping Proportions.

		Low skipping proportion		High skipping proportion
Correlation		NMAR	MAR	NMAR	MAR
$\mathrm{\rho }=0$	DIC	20,278.23	20,277.16	24,192.70	24,194.22
	LPML	−10,625.71	−10,622.35	−12,639.62	−12,640.24
$\mathrm{\rho }=-0$ .5	DIC	20,077.83	20,373.61	24,100.94	24,563.02
	LPML	−10,481.75	−10,582.50	−12,551.96	−12,680.42
$\mathrm{\rho }=-0$ .8	DIC	20,245.78	20,577.46	24,137.64	24,578.39
	LPML	−10,563.83	−10,677.20	−12,566.58	−12,722.20

In conclusion, the results demonstrate that under various conditions, the NMAR model provides a better fit to the data compared to the MAR model, as evidenced by DIC and LPML values.

Empirical Example

We analyzed the data from the computer-based PISA 2018 mathematics cognitive test in Albania. According to the 2018 PISA mathematics assessment framework, four attributes belonged to mathematical content knowledge: change and relationship $\left({\mathrm{\alpha }}_1\right)$ , quantity $\left({\mathrm{\alpha }}_2\right)$ , space and shape $\left({\mathrm{\alpha }}_3\right)$ , and uncertainty and data $\left({\mathrm{\alpha }}_4\right)$ . The item IDs and Q matrices are shown in Table 10; these items were dichotomous. Item responses were coded 0 (no credit), 1 (full credit), 6 (not reached), 7 (not applicable), 8 (invalid), or 9 (nonresponse). Because the nonresponses were discrete, we treated code 9 as skipped/omitted responses. Furthermore, with some examinees’ responses removed due to being coded as not applicable or invalid, 798 examinees remained in the data set. Furthermore, another 40 examinees coded 6 (not-reached) were removed. The final sample size is 758. The overall skipping proportion was 4.9%, and the skipping proportions of each item ranged from 0% to 14.8%.

Table 10.

The Q Matrix in the Real Data.

Attribute	CM033Q01	CM474Q01	CM155Q01	CM155Q04	CM411Q01	CM411Q02	CM803Q01	CM442Q02	CM034Q01
${\mathrm{\alpha }}_1$	0	0	1	1	0	0	0	0	0
${\mathrm{\alpha }}_2$	1	0	0	0	0	0	0	0	1
${\mathrm{\alpha }}_3$	0	1	0	0	1	0	0	1	0
${\mathrm{\alpha }}_4$	0	0	0	0	0	1	1	0	0

Tables 11 and 12 show the parameter estimates and standard deviations of the person, item, and attribute parameters. The results indicate that the estimated correlation coefficient between the higher-order ability and skipping ability parameter ${\sigma }_{{\theta }^h{\theta }^s}$ was negative. This suggests that as the higher-order ability decreased, the latent person parameter associated with skipping increased, indicating that examinees with higher higher-order abilities found it more challenging to skip items. The positive estimates for the attribute slope parameters imply that as the higher-order ability parameter ${\mathit{\theta }}_{\mathit{i}}^{\mathit{h}}$ increases, examinees demonstrate a better mastery of the knowledge attribute, which aligns with real-world expectations. Furthermore, the estimated item mean parameter ${\mu }_{\beta }$ was −2.067, which implies that the mean guessing probability was approximately .11. The correlation coefficient between the item intercept parameter ${\beta }_j$ and interaction parameter ${\delta }_j$ was negative ( ${\sigma }_{\beta \delta }<0$ ). In addition, all item intercept parameters ${\beta }_j$ were negative, which suggests that the guessing probability of all the items was lower than .5. All item interaction parameters ${\delta }_j$ were positive, which satisfied $g_j<1-s_j$ , as expected. The smaller the skipping difficulty parameter ${\beta }_j^s$ , the higher the skipping proportion for this item; for example, for CM155Q01, the estimation of ${\beta }_j^s$ was 2.621 and the skipping proportion was 14.8%. For CM442Q02, the estimation of ${\beta }_j^s$ was 2.998 and the skipping proportion was 11.6%. Although the overall skipping proportion in the real data was 4.9%, certain items exhibited notably high skipping proportions, indicating the need to consider them as missing data.

Table 11.

Estimates and Standard Deviations of the Parameters for the Real Data.

Statistics	${\mu }_{\beta }$	${\mu }_{\delta }$	${\mu }_{\beta }^s$	${\sigma }_{\beta }^2$	${\sigma }_{\delta }^2$	${\sigma }_{{\beta }^s}^2$	${\sigma }_{\beta \delta }$	${\sigma }_{\beta {\beta }^s}$	${\sigma }_{\delta {\beta }^s}$	${\sigma }_{{\theta }^h{\theta }^s}$
Est.	−2.067	2.471	4.787	2.543	1.537	3.302	−0.870	0.475	−0.634	−1.122
SD	0.475	0.405	0.560	1.688	1.208	2.093	1.141	1.206	0.997	0.142
Statistics	${\sigma }_{{\theta }^s}^2$	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\lambda }_4$	${\gamma }_1$	${\gamma }_2$	${\gamma }_3$	${\gamma }_4$
Est.	3.076	−0.218	0.450	0.779	1.414	2.416	2.396	2.171	1.702
SD	0.568	0.397	0.525	0.398	0.653	0.572	0.609	0.416	0.357

Note. Est. = estimated value; SD = standard deviation of the posterior distribution.

Table 12.

Estimates and Standard Errors of the Item Parameters for the Real Data.

Parameter	Statistics	033Q01	474Q01	155Q01	155Q04	411Q01	411Q02	803Q01	442Q02	034Q01
${\beta }_j$	Est.	−0.157	−0.532	−1.546	−1.452	−1.837	−1.354	−3.913	−4.430	−3.386
	SD	0.246	0.145	0.537	0.198	0.270	0.182	0.550	0.535	0.734
${\delta }_j$	Est.	2.440	1.618	4.219	1.344	2.573	0.825	3.203	2.990	2.801
	$\mathrm{SD}$	0.499	0.250	0.829	0.257	0.375	0.369	0.632	0.588	0.797
${\beta }_j^s$	Est.	5.907	6.807	2.621	4.637	3.665	5.026	4.887	2.998	7.808
	SD	0.402	0.575	0.185	0.289	0.233	0.325	0.311	0.199	0.904

Note. Est = estimated value; SD = standard deviation of the posterior distribution.

Discussion

This study addresses the prevalent issue of item-level missing data in cognitive diagnostic assessments, where examinees skip certain items. Such missing data poses a risk of providing inaccurate diagnostic feedback and may lead to biased outcomes when compared to those with complete data. Therefore, this study aims to develop a missing data model for skipped items, ultimately providing more accurate cognitive feedback for examinees. This not only facilitates a deeper comprehension of individual cognitive processes but also furnishes more robust data support for educational policies and practices.

Our proposed model combines the DINA model, which captures responses, and the 1PLM, which represents the missing indicator model. A bivariate normal distribution for person parameters and a multivariate normal distribution for item parameters were used to capture the relationships between these two models in a higher-level framework. Our model not only strives for more precise knowledge structures derived from diagnostic feedback but also sheds light on examinees’ response behaviors, offering valuable insights into their cognitive processes.

We employed the Bayesian MCMC algorithm to estimate the model parameters. Through comprehensive simulation studies, we found our model consistently yielded more accurate parameter estimates compared to other methods. These simulation studies also facilitated an in-depth evaluation of the model’s performance, considering various missing mechanisms, proportions of skipped items, and correlations between higher-order ability and potential skipping tendency. Under these diverse simulation scenarios, our model consistently demonstrated reliable estimation accuracy. Furthermore, the HO-DINA model, which characterizes examinees’ cognitive attributes, further enhanced the accuracy of our parameter estimates. Notably, we observed substantial improvements in ACCRs and PCCRs, specifically, ACCRs and PCCRs were significantly higher in the NMAR model, indicating that our model better fits non-ignorable missing data mechanisms. To demonstrate the practical applicability of our model, we successfully applied it to the PISA 2018 computer-based test.

While this study offers valuable insights, it comes with certain limitations, suggesting potential directions for future research. Firstly, it’s essential to consider time constraints as a factor contributing to item skipping. Secondly, we have not explored other potential reasons for skipped items in cognitive assessments. Future research should delve into these additional factors. Lastly, cognitive assessments frequently exhibit various types of item nonresponses, including not-reached and skipped items. Future research could refine our model to address these diverse missing data scenarios.

Conclusion

In this study, we introduced a novel model to address skipped items in cognitive diagnostic assessments, combining the missing data model with the item response model under a cohesive hierarchical structural framework. By leveraging the Bayesian MCMC algorithm, robust parameter estimation results are obtained.

Extensive simulations demonstrated the superior performance of our model across a range of scenarios, consistently producing highly accurate parameter estimates. Notably, our model exhibited a significant advantage in ACCRs and PCCRs when compared to the NMAR and MAR models, particularly in the presence of non-ignorable missing data.

The validation of our model using the 2018 PISA mathematics dataset further highlights its practical effectiveness. Our research has profound implications for cognitive diagnostic assessments, offering a transformative approach that addresses skipped items comprehensively. It lays the foundation for potential advancements in the field, promising more accurate and insightful assessments in the future.

Supplemental Material

sj-docx-1-sgo-10.1177_21582440241288005 – Supplemental material for Modeling Skipped Items in Cognitive Diagnostic Assessments

Supplemental material, sj-docx-1-sgo-10.1177_21582440241288005 for Modeling Skipped Items in Cognitive Diagnostic Assessments by Lidan Liang, Jing Lu, Jiwei Zhang and Guanghui Cheng in SAGE Open