Glas and Pimentel (2008) proposed a sequential model with a linear restriction on the item parameters to model the not-reached items. Specifically, the missing indicator matrix D with element d_ij is given by:

where d_ij = 1indicates that the ith examinee drops out the jth item. Because of the small overall proportion of not-reached responses, the appropriate model must have few parameters to be estimable (Lord, 1983). The one-parameter logistic model (1PLM; Rasch, 1960) is adopted to model the missing indicators, thus the dropping-out probability of examinee i on item j is:

and

where βjd represents the so-called item difficulty parameter for item j, andθid denotes the ith examinee’s dropping-out propensity. Also, βjd=η0 when j = J, where η₀ is the difficulty threshold of the last item, and η₁ models a uniform change in the probability as a function of the item position in the test. Usually, the parameter η₁ is negative, and hence it is more likely to drop out the test at later position items of the test.

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

where s_j and g_j are the slipping and guessing probabilities of the jth item, respectively, 1−s_j−g_j = IDI_j is the jth item discrimination index (de la Torre, 2008), and α_ik is the kth attribute of the ith examinee, with α_ik = 1 if examinee i masters attribute k and α_ik = 0 if examinee does not master attribute k. The Q matrix (Tatsuoka, 1983) is an J×K matrix, with q_jk,q_jk = 1 denoting that the attribute k is required for answering the jth item correctly and q_jk = 0 if the attribute k is not required for answering the jth item correctly.

Equation (4) can be reparameterized as the reparameterized DINA model (DeCarlo, 2011).

In addition, logit(x)=log(x1-x), thus Equation (4) can be reformed as,

where β_j and δ_j are the item intercept and interaction parameter, respectively, and they are assumed to follow a bivariate normal distribution as follows:

The higher-order structure is very flexible because it can reduce the number of model parameters and can provide higher-order abilities and more accurate attribute structures. Because the attributes in a test are often correlated, the higher-order structure (de la Torre and Douglas, 2004; Zhan et al., 2018) for the attributes is expressed as,

where P(α_ik = 1) is the probability that the ith examinee masters the kth attribute, θih is the higher-order ability of examinee i, and γ_k and λ_k are the slope and intercept parameters of attribute k, respectively. The slope parameter γ_k is positive because the knowledge attribute is mastered better with the increased ability θih.

If the observation probability p(y_ij|d_ij, β_j, δ_j, α_ik) does not depend on θ^d, when θ^h and θ^d are independent, then the missing data are ignorable. In this situation, this model is treated as a MAR model. Let p(y_ij|d_ij, β_j, δ_j, α_ik) be the measurement model for the observed data. In addition, let p(dij|θid,η0,η1) be the measurement model for the missing data indicators, and p(θ^h) and p(θ^d) are densities of θ^h and θ^d, respectively. To model non-ignorable missing data, it is assumed that θih and θidfollow a bivariate normal distribution N(μ_P, Σ_P); thus, the two models describe the two missing mechanisms (i.e., MAR and NMAR). Next, we introduce the two missing data models for the not-reached items.

In Equations (2) and (9), the linear parts of 1PLM and the HO-DINA model can be written as follows:

To eliminate the trade-offs between ability θid and dropping-out threshold parameter βjd and between the higher-order ability person parameter θ^h and the attribute intercept λ_k, the mean population level of person parameters is set to zero, that is, μθh=0andμθd=0. σθh= 1 is fixed to eliminate the scale trade-off between θih and γ_k (Lord and Novick, 1968; Fox, 2010). In addition to the identifications, two local independence assumptions are made, that is, the α_ik values are conditionally independent given θih, and the Y_ij values are conditionally independent given α_i.

In the Bayesian framework, 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 and Eddy, 1979; Ibrahim et al., 2001) are used to compare the differences in the missing mechanism models according to the results of MCMC sampling. Let,

The DIC is given by,

On the basis of the posterior distribution of Dev(Y,D,Ω), the DIC was defined as,

where Dev¯=E(Dev(Y,D,Ω)|Y,D)≅1R∑r=1R[cpsbreak]Dev(Y,D,Ωr), which is the posterior mean deviance and is a Bayesian measure of fit, r = 1,…,R denotes the rth iteration of the algorithm, and Dev^=Dev(Y,D,Ω¯), which is the effective number of parameters, is a Bayesian measure of complexity, with Ω¯=E(Ω|Y,D)≅1R∑r=1RΩr. A smaller DIC indicates a better model fit.

The conditional predictive ordinate (CPO) index of the two models was computed. Let Q_ij,max = max_1≤r≤R{−logf(Y_ij,D_ij|Ω^r)}. Thus,

The summary statistic for log(CPOij^) is the sum of their logarithms, which is termed the LPML and is given by,

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

In the three simulation studies, the item parameters were sampled from the following distributions: (β_jδ_j)∼MVN((μ_βμ_δ), Σ_I), μ_β = −2.197,μ_δ = 4.394, Σ_I = (1−0.8−0.81). These values were used in Shan and Wang (2020) study. The dropping-out proportions across three levels (i.e., low, medium, and high) were varied by setting different combinations of η₀ and η₁. That is, the dropping-out proportion was 3.8 (low) when η₀ = 1,η₁ = −0.7; the dropping-out proportion was 12 (medium) when η₀ = 1,η₁ = −0.32; and the dropping-out proportion was 25% (high) when η₀ = 1,η₁ = −0.18.

The attribute intercept parameters were λ=(−1,−0.5, 0, 0.5, 1), and the attribute slope parameters wereγ_k = 1.5 for all attributes, which were consistent with those in the study by Shan and Wang (2020). Three Q matrices with different numbers of attributes (Figure 1) were considered, and the three Q matrices were taken from Xu and Shang (2018) study and Shan and Wang (2020) study.

The person parameters θih and θid were simulated from the bivariate normal distribution (θhθd)∼MVN((00),(1σθhθdσθhθdσθd2)),whereσθd2= 0.25. Three levels of correlation between θih and θid were considered for ρθihθid:0 (uncorrelated), −0.5 (medium), and −0.8 (high). The missing data due to dropping-out items were simulated in the following manner. The three levels of dropping-out proportions were 3.8% (low), 12% (medium), and 25% (high).

The priors of η₀ and η₁ were η₀ ∼ N(0,2) and η₁ ∼ N(0,2), respectively. The priors of the item parameters β_j and δ_j were assumed to have a bivariate normal distribution: (β_jδ_j)∼N((μ_βμ_δ), Σ_I). The priors of the person parameters were assumed to follow a bivariate normal distribution: (θ^hθ^d)∼N((00), Σ_P). The priors of the higher-order structure parameters were expressed as λ_k∼N(0,4) and γ_k ∼ N(0,4)I(γ_k > 0), the priors of the covariance matrix of the person were expressed as σ_θ^hθ^d∼U(−1,1) and σθd2∼Inv-(2,2), the priors of the covariance matrix of the item parameters were expressed as Σ_I Inv-Wishart (ΣI0-1,vI0), and the hyperpriors were specified as Σ_I0 = (1001), v_I0 = 2,k_I0 = 1,μ_β∼N(−2.197,2), and μ_δ∼N(4.394,2)I(μ_δ > 0). The hyperpriors specified above were on a logit scale for β and δ and were consistent with those reported by Zhan et al. (2018). The mean guessing effect was set at 0.1, which was roughly equal to a logit value −2.197 for μ_β. A standard deviation of 2 on the logit scale for μ_βindicated that the simulated mean guessing effect changed from 0.026 to 0.314. In addition, the mean slipping effect was also set at 0.1, which indicated that μ_δ was approximately 4.394 on the logit scale. The simulated mean slipping effect changed from 0.007 to 0.653 under a standard deviation of 2 on the logit scale for δ.

The initial values of the model parameters were as follows: β_j = 0, δ_j = 0for j = 1,…,J,θih=0, θid=0for i = 1,…,N, σ_θ^hθ^d = 0, σθd2=1, η₀ = 0,η₁ = 0,μ_β = 0,μ_δ = 0, Σ_P = (1001), μ_P = (00), and σθh2=1. In addition, λ_k = 0,γ_k = 1for k = 1,,K, and α=(α₁₁⋯α_1K⋮⋱⋮α_N1⋯α_NK), where α_ik (i = 1,…,N,k = 1,…,K) were sampled from 0 to 1 randomly. The proposal variances were chosen to give Metropolis acceptance rates between 25% and 40%. The Markov chain length was set at 10,000 so that the potential scale reduction factor (PSRF; Brooks and Gelman, 1998) was less than 1.1 for all parameters, which implied proper chain convergence. Five thousand iterations were treated as burn-in. The final parameter estimates were obtained as the average of the post-burn-in iterations.

In terms of evaluation criteria, the bias and root mean squared error (RMSE) are used to assess the accuracy of the parameter estimates. In particular, the bias for parameter η was,

and the RMSE for parameter η is defined as,

where η is the true value of the parameter, and η^(r) is the estimate for the rth replication. There were R = 30 replications for each simulation condition. The recoveries of attributes are evaluated using the attribute correct classification rate (ACCR) and the pattern correct classification rate (PCCR):

where I(αik^=αik) is the indicator function that is, I(αik^=αik)=1 if αik^=αik, otherwise I(αik^=αik)=0.

In simulation study I, the different numbers of examinees and items were considered to estimate the model parameters under a fixed number of five attributes. Three conditions were considered in this simulation: (a) 500 examinees and 30 items, (b) 1,000 examinees and 30 items, and (c) 500 examinees and 20 items. The correlation between θih and θid was −0.3, and the dropping-out proportion was medium.

Table 1 presents the bias and RMSE of the ability parameters and item parameters, as well as the attribute parameter estimates. For the 30 items and the 5 attributes (please see the first four columns of Table 1), the item parameter estimates improve when the number of examinees increases from 500 to 1,000, the bias and RMSE of δ and μ_βdecrease, and the RMSE of β,μ_δ, and item covariance matrix elements reduce. For the 500 examinees and the 5 attributes (please see the middle four columns of Table 1), the person parameter estimates improve when the number of items increases from 20 to 30, and θ^h and σθd2 are more accurate. The ACCRs and PCCRs are presented in Table 2. The ACCRs and PCCRs could be recovered satisfactorily with a larger sample and longer test length. The ACCRs and PCCRs decrease when the number of examinees or test length decreases (please see the first three columns in Table 2), and the changes are particularly marked when the test length is reduced. Figure 2 shows the PSRF of several items and attribute parameters under 500 examinees and 30 items. It is observed that the item intercept parameter β, the interaction parameter δ, the attribute slope parameter γ, and the attribute intercept parameter λ converge at 5,000 iterations, and the convergence of β and δ are significantly faster than that of λ and γ.

This simulation study was conducted to investigate the parameter recovery of different numbers of attributes for fixed 500 examinees and 20 items. The correlation between θih and θid was set at −0.3, and the dropping-out proportion was medium.

The last four columns of Table 1 show the results of simulation study II. The RMSE of the estimates of item and person parameters with attribute K = 5 are smaller than those with attribute K = 3. The RMSE of the attribute slope parameters and intercept parameters recover more satisfactorily with attribute K = 3 than with attribute K = 5. The last two columns of Table 2 show the ACCRs and PCCRs for simulation study II. The ACCRs with attribute K = 3 are higher than those with attribute K = 5 and improve from 0.957 to 0.987 on average. Moreover, the PCCRs are significantly higher when the number of attributes decreases. That is, the PCCR with attribute K = 5 is 0.811, and the PCCR with attribute K = 3 is 0.961.

The purpose of this simulation study was to investigate the parameter recovery with the NMAR model, MAR model, and HO-DINA model that ignores the not-reached items under different simulation conditions. The data were generated using the proposed model with the NMAR mechanism. A total of 500 examinees answered 30 items, and each item had 5 attributes. Three dropping-out proportions (i.e., 3.8% [low], 12% [medium], and 25% [high]) and three correlations between θih and θid (i.e., 0 [uncorrelated], −0.5 [medium], and −0.8 [high]) were manipulated. Thus, there were 3 × 3 simulation conditions.

Table 3 shows the bias and RMSE of the parameters of three models with low dropping-out proportions under different correlations between θih and θid. Results show that the parameter estimates from the three models are similar when the correlation between θih and θid is 0. When the correlation between θih and θid increases, the bias and RMSE of η₁, β, Σ_I, and γ in the NMAR model are much smaller than those in the MAR and HO-DINA models. Moreover, for low dropping-out proportions, when the correlation between θih and θid increases, the bias of the person parameters of the three models changes very little, whereas the RMSE of the person parameters in the MAR and HO-DINA models increases significantly. As expected, the NMAR model has higher accuracy of parameters than that of the other two models. Furthermore, the parameter estimates of the MAR and HO-DINA models are similar for all simulation conditions because θih and θid are uncorrelated in both the MAR and HO-DINA models, which ignore the not-reached items. Table 4 shows the bias and RMSE of the parameters of the three models with medium dropping-out proportions under different correlations between θih and θid. Similar parameter estimates are obtained from the three models when the correlation between θih and θid is 0. When the correlation between θih and θid increases, not only the bias but also the RMSE of the person parameters are lower in the NMAR model than those in the MAR and HO-DINA models, and the other results are similar to those with low dropping-out proportions. Table 5 shows the bias and RMSE of the parameters of the three models with high dropping-out proportions under different correlations between θih and θid. We find that the parameter estimates improve significantly with high dropping-out proportions. Figure 3 shows the bias of the estimates of item mean vector and the item covariance matrix elements in the NMAR and MAR models under different dropping-out proportions and correlations between θih and θid. The results show that the estimates of the parameters are more accurate in the NMAR model than those in the MAR model when the correlation is increased. Moreover, it is observed that the bias of the parameters of the NMAR model is close to 0 as the correlation between θih and θid increases. In contrast, the bias of the parameters of the MAR model is significantly larger than that of the NMAR model. Figure 4 shows the RMSE of the estimates of the item mean vector and the item covariance matrix elements in the NMAR and MAR models under different dropping-out proportions and correlations between θih and θid. The results show that the RMSE of the item mean vector in the NMAR model improves slightly than that in the MAR model. Moreover, the RMSE of the item covariance matrix elements shows significant improvements, and the estimates of the item covariance matrix elements are precise when the correlation is high. Figure 5 shows the ACCRs and PCCRs under nine simulation conditions. Detailed results are provided in Supplementary Table 1. It is found that ACCRs and PCCRs in the NMAR model are improved significantly when the missing proportion or the correlation between θih and θid is high. This indicates that the MAR model could not recover the attribute pattern effectively when the missing data mechanism is indeed non-ignorable. Table 6 shows the model selection results. The differences in DIC and LPML are not obvious when the correlation between θih and θid is 0. The DICs of the NMAR model are smaller than those of the MAR model under nine simulation conditions. Moreover, the LPMLs of the NMAR model are higher than those of the MAR model. Thus, the DIC and LPML indices are able to select the true model accurately.