GDINA : An R Package for Cognitive Diagnosis Modeling

Cognitive diagnosis models (CDMs) have attracted increasing attention in educational measurement because of their potential to provide diagnostic feedback about students’ strengths and weaknesses. This article introduces the feature-rich R package GDINA for conducting a variety of CDM analyses. Built upon a general model framework, a number of CDMs can be calibrated using the GDINA package. Functions are also available for evaluating model-data ﬁt, detecting diﬀerential item functioning, validating the item and attribute association, and examining classiﬁcation accuracy. A grapical user interface is also provided for researchers who are less familar with R . This paper contains both technical details about model estimation and illustrations about how to use the package for data analysis. The GDINA package is also used to replicate published results, showing that it could provide comparable model parameter estimation.


Introduction
Cognitively diagnostic assessments (CDA; de la Torre and Minchen 2014; Nichols, Chipman, and Brennan 1995) have gained increasing popularity in the past decades in the field of educational measurement. Traditional standardized educational assessments usually base on unidimensional item response theory (IRT; e.g., de Ayala 2013), which assumes that a single latent trait (or overall ability) is measured. Consequently, students are located on a continuum based on their performance in assessments using appropriate IRT models. In contrast, CDAs usually depend on cognitive diagnosis models (CDMs) with an intention to provide diagnostic information about students' strengths and weaknesses. To obtain such information, the assessments are typically designed to measure a set of finer-grained skills, which are usually referred to as attributes, and treated as binary latent variables with outcome 1 for mastery and 0 for nonmastery. The major goal of CDM analyses is to infer students' attribute profiles from their item responses.
A wide array of CDMs (e.g., Rupp, Templin, and Henson 2010) have been developed, and most of them, if not all, consist of two major components. One component is an item and attribute association matrix, or Q-matrix (Tatsuoka 1983), which specifies whether an attribute is needed to perform an item correctly. The Q-matrix is usually developed by domain experts, and assumed to be correct in the following CDM analyses. The second component is referred to as the condensation rule (Maris 1999), which defines how attributes are "condensed" to yield manifest item responses. Some CDMs are formulated to accommodate some specific condensation rule. For example, the deterministic inputs, noisy "and" gate (DINA; Haertel 1989) model employs a conjunctive rule, and assigns the highest probability of answering correctly to individuals that possess all of the required attributes. The deterministic inputs, noisy "or" gate (DINO; Templin and Henson 2006) model, however, adopts a disjunctive rule, and assigns the highest probability of answering correctly to individuals mastering at least one of the required attributes. In addition, the additive CDM (A-CDM; de la Torre 2011) assumes that each required attribute contributes to the success probability uniquely and independently. Apart from these specific CDMs, general or saturated CDMs subsuming many widely-used specific CDMs have also been developed, including the generalized DINA (G-DINA; de la Torre 2011) model, the general diagnostic model (GDM;von Davier 2008) and the log-linear CDM (LCDM; Henson, Templin, and Willse 2009).
The present work introduces the R (R Core Team 2020) package GDINA (Ma and de la Torre 2020b) for conducting a variety of CDM analyses based on the G-DINA model (de la Torre 2011) and its extensions. The package is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=GDINA. Currently, there exist a few packages for the R programming environment which implement CDM analyses, including the ACTCD (Chiu and Ma 2018), CDM (George, Robitzsch, Kiefer, Gross, and Ünlü 2016), dina (Culpepper 2015), and NPCD (Zheng and Chiu 2019) packages. The dina package estimates the DINA model using the Gibbs sampler, and the ACTCD and NPCD packages are designed for nonparametric CDM analyses. The CDM package is capable of estimating the GDM, the structured CDMs (Formann 1992), the regularized latent class model (Chen, Li, Liu, and Ying 2017), as well as the G-DINA model. In addition to these R packages, some general purpose commercial software programs, such as Mplus (Muthén and Muthén 2017) and Latent GOLD (Vermunt and Magidson 2016), are also capable of fitting the G-DINA model by imposing some specific constraints. However, since they are not designed for CDM analyses, they usually lack some important functionalities (e.g., Q-matrix validation), and tend to be slow.
The GDINA package is primarily developed for conducting CDM analyses using the G-DINA model (de la Torre 2011) and its extensions (e.g., Chen and de la Torre 2013; Ma and de la Torre 2016; Ma, Terzi, Lee, and de la Torre 2017; Ma 2019b), hence the name. It has several distinguishing features compared with existing R packages for cognitive diagnosis at the time of writing. In particular, the GDINA package provides a set of unique tools for model diagnostics under the G-DINA model framework. For example, the GDINA package is the only program that allows for empirically validating the Q-matrix under a general model framework using a variety of approaches (de la Torre and Chiu 2016; de la Torre and Ma 2016; Ma and de la Torre 2020a; Najera, Sorrel, and Abad 2019) and selecting models at item level using the Wald test, likelihood ratio test or Lagrange multiplier test. The GDINA package also allows users to assess global model-data fit using the M 2 statistic for dichotomous response (Hansen, Cai, Monroe, and Li 2016;Liu, Tian, and Xin 2016) and the M ord statistic for ordinal response (Ma 2020). In addition, the GDINA package allows users to calibrate the generalized multiple-strategy models for dichotomous response (Ma and Guo 2019) and diagnostic tree model for polytomous response (Ma 2019b) when multiple strategies exist. Furthermore, the GDINA package provides a routine to validate the Q-matrix, choose the most appropriate CDMs for each item, and calibrate the selected CDMs at one fell swoop. Last, the GDINA package is the only R package that offers a graphical user interface for various CDM analyses.
The paper is organized as follows. In Section 2 we review the G-DINA model framework, including various parameterizations of item response functions and joint attribute distribution. In Section 3, we present the details of the EM algorithm for structural parameter estimation. We introduce the features of the GDINA package in Section 4 and analyze a real dataset for illustrative purpose in Section 5. In Section 6, we compare the GDINA package and commercial software programs in terms of parameter estimation and speed based on two sets of real data. We conclude with a discussion on possible future developments in Section 7.

The G-DINA model framework
Suppose a test with J items measures K binary attributes. The association between items and attributes is specified in a J × K Q-matrix (Tatsuoka 1983), with element q jk = 1 indicating item j measures attribute k and q jk = 0 indicating item j does not measure attribute k. Also, K binary attributes produce 2 K attribute profiles, each labeling a unique latent class. The attribute profile for latent class c is denoted as α c = [α c1 , . . . , α cK ] , where α ck = 1 if attribute k is mastered and α ck = 0 if not. Superscript is used to denote transposition. Let a = [a 1 , . . . , a i , . . . , a N ] be a vector of attribute profiles of N individuals in a sample, where a i ∈ {α 1 , . . . , α c , . . . , α 2 K } is a K-dimensional random vector representing the attribute pattern of individual i. Let Y ij be a response variable of individual i to item j, following a Bernoulli distribution with an observed realization y ij . Also, denote Y i = {Y ij } as item response vector of individual i and Y = [Y 1 , . . . , Y N ] the item response matrix from N individuals. To formulate a CDM, we need to define the item response function for the relation between item response and attribute profiles, and the model for joint attribute distribution, which gives the population proportion of each latent class.

Item response function
When not all attributes are required for item j, the G-DINA model (de la Torre 2011) collapses 2 K latent classes into 2 K * j latent groups, where K * j = K k=1 q jk . The collapsed latent classes have the same success probability to item j. To simplify the notation, the first K * j attributes are assumed to be required for item j, and α * lj is the reduced attribute vector consisting of the columns of the required attributes, where l = 1, . . . , 2 K * j . The probability of individual i with attribute profile α c answering item j correctly is denoted by P(Y ij = 1|α c ) = P j (α c ). Note that if latent class c is collapsed into latent group l, we have P j (α c ) = P(α * lj ), where for the latter, the subscript j is dropped to avoid the redundancy. The item response function (IRF) of the G-DINA model is given by where g ·] stands for identity, log or logit link functions. Additionally, δ j0 is the intercept for item j, δ jk is the main effect due to α k , δ jkk is the interaction effect due to α k and α k , δ j12···K * j is the interaction effect due to α 1 , . . ., α K * j . By setting appropriate constraints, as shown in de la Torre (2011), several widely used CDMs can be obtained. Specifically, to obtain the DINA model, all terms in the identity link G-DINA model, except δ 0 and δ 12...K * j , are constrained to zero, that is, For the DINO model, the IRF is given by This results in two parameters for the DINO model as well as the DINA model for each item, regardless of the number of attributes required. In addition, by setting all interactions at zero in the G-DINA model, the A-CDM, linear logistic model (LLM; Maris 1999), and reduced reparameterized unified model (R-RUM; Hartz 2002) can be obtained. Specifically, the A-CDM is the constrained identity-link G-DINA model without the interaction terms. It can be written by The LLM is the logit link G-DINA model without any interaction terms, and can be formulated as The R-RUM is the log-link G-DINA model without any interaction terms, and can be given by These three models assume that mastering attribute α lk raises the probability of success on item j, but the increases may be different due to the different scales they employ. As additive models, mastery of each attribute does not affect other attributes. Additionally, they all have the same number of parameters for item j (i.e., K * j + 1). Maris (1999) discussed a variety of approaches to parameterizing the joint attribute distribution. In this paper, the probability mass function of attribute profile α c is denoted by

Joint attribute distribution
where λ is a vector of unknown structural parameters. One of the simplest approaches is the independent model, which assumes that all attributes are independent statistically and can be written by The independent model involves K parameters, λ = [P(α 1 = 1), . . . , P(α K = 1)] , where P(α k = 1) is the probability of mastering attribute k. Another model with simple parameterization is the saturated model, which involves 2 K parameters, λ = [π 1 , . . . , π 2 K ] , with the constraint of 2 K c=1 π c = 1. Note that π c is sometimes referred to as mixing proportion parameter, indicating the proportion of individuals in latent class c. The saturated model is considered the most general parameterization of the joint attribute distribution, but tends to involve too many parameters when K is large. In addition, Xu and von Davier (2008) proposed to use a loglinear model to smooth the joint attribute distribution, which can be written as where n(α c ) = N π c is the number of individuals with attribute profile α c . The above loglinear structural model considers main effects and first-order interactions. It has 1 + K(K + 1)/2 parameters, namely, λ = [λ 0 , λ 1 , . . . , λ K−1,K ] . It is very flexible in that it can be simplified by removing all interactions between attributes, or extended by including higher-order attribute interactions, but its parameters do not have straightforward interpretation. Another method for parameterizing the joint attribute distribution is the higher-order model (de la Torre and Douglas 2004), which assumes that the mastery of each attribute is influenced by a higher-order latent trait and that their relation is defined using IRT models, such as the two parameter logistic (2PL; see de Ayala 2013) model: where θ represents the unidimensional general ability and λ = [λ 01 , . . . , λ 0K , λ 11 , . . . , λ 1K ] is a vector of higher-order structural parameters. Note that setting λ 1k = λ 1 ∀k yields the one parameter logistic (1PL; see de Ayala 2013) model, and setting λ 1k = 1 ∀k produces the Rasch model (Rasch 1960). Under the assumption of local independence, and the probability mass function can be obtained by integrating out θ, which can be approximated by Gauss-Hermite quadrature: whereθ s and W (θ s ) are the quadrature nodes and weights, respectively. The accuracy of the approximation can be improved by increasing the number of quadrature nodes S.
Note that the higher-order model defined in de la Torre and Douglas (2004) is more general by allowing the higher-order latent trait to be multidimensional, but only the unidimensional θ is considered in this paper, as well as in the GDINA package. There is no existing software program that can accommodate higher-order CDMs with multidimensional higher-order latent traits so far. This is one of the features that we consider incorporating into the GDINA package later.

Model estimation
Let γ = [δ , λ ] denote a vector of structural parameters involved in a CDM, including item parameters δ and joint attribute distribution parameters λ. In addition, for individual i, attribute profile a i , as well as higher-order ability θ i when a higher-order model is employed for joint attribute distribution, need to be estimated, which are typically referred to as incidental parameters. The G-DINA model was estimated by maximizing the marginalized likelihood (de la Torre 2011), which, as shown in Bock and Aitkin (1981), can be implemented using the expectation-maximization algorithm (EM; Dempster, Laird, and Rubin 1977). The EM algorithm is a general procedure for finding MLEs when missing data exist. In the current context, for individual i, item responses y i can be regarded as "incomplete" data, person parameter a i as missing data, and x i = [y i , a i ] as the "complete" data. When a higherorder model is used, the complete-data for individual i is The goal of the EM algorithm is to find the maxima of the incomplete-data likelihood indirectly by maximizing the complete-data log likelihood iteratively. More specifically, the EM algorithm consists of two steps: the expectation (E) step and the maximization (M) step. In the E-step, we need to calculate the so-called Q function (Dempster et al. 1977), which is the expected log likelihood of the complete-data conditional on the observed data and current parameter estimates, and takes the following form: where γ denotes the parameter estimates from the previous step. In the M-step, the Q function is maximized. The E-and M-steps repeat until certain convergence criteria have been met.
It can be shown that, when individuals are independent and the joint attribute distribution is not modeled using the higher-order model, where n c is the expected number of individuals in latent class c and r jc is the expected number of individuals in latent class c who answer item j correctly. They can be calculated by where P(α c |y i , γ ) is the posterior probability of individual i being assigned to latent class c, and can be calculated using the Bayes rule: Here, When individuals are independent and the joint attribute distribution is parameterized using a higher-order model, where n s is the expected number of individuals having abilityθ s , which can be calculated by and r ks is the expected number of individuals with abilityθ s mastering attribute k, and can be calculated by In the E-step, n c and r jc , as well as n s and r ks when a higher-order model is employed, are calculated based on γ , and then, in the M-step, Q(γ; γ ) is maximized with respect to model parameters γ.

Implementation in the GDINA package
The main function of the R package GDINA is GDINA(), which allows the calibration of a variety of CDMs within the G-DINA model framework. In the GDINA package, the G-DINA model is specified similar to the generalized linear model using design matrix and link function (de la Torre 2011). In particular, let δ j = [δ j0 , . . . , δ j12...K * j ] be a vector of parameters of item j and P j = {P(α * lj )} be a vector of item success probabilities. The G-DINA model is written as where M j is the design matrix (de la Torre 2011) and g[·] the link function. Let us assume item j requires 2 attributes, for the G-DINA model, Note that the aforementioned DINA, DINO, A-CDM, LLM and R-RUM models can also be written in this form by specifying appropriate design matrices. For example, the M j for the DINA and DINO models can be written by respectively. Furthermore, it is possible to define the design matrix and link function to obtain new models. For example, as pointed out by de la Torre (2011) From the design matrices, it is straightforward to note that the Bug-DINA and Bug-DINO models are just reparameterized DINA and DINO models. In the GDINA package, one can define their own models by specifying the design matrix and link function for each item.
Model parameters are estimated using the EM algorithm. The E-step and many other computationally intensive functions are written in C++ through the Rcpp (Eddelbuettel and François 2011) and RcppArmadillo (Eddelbuettel and Sanderson 2014) R packages to speed up the execution. In the M-step, for a few models, closed-form solutions exist; whereas for other models, some optimization routines are needed to maximize the Q function. Possible solvers include optim from the stats package (R Core Team 2020), slsqp fromthe nloptr package (Johnson 2010), auglag from the alabama package (Varadhan 2015), or solnp from the Rsolnp package (Ghalanos and Theussl 2015). By default, the solvers are automatically chosen according to the model to be calibrated, and if one optimization routine fails, another one may be employed. Users can specify the solver in the GDINA() function if they have any preference.
Several methods are available to further analyze the object returned from the GDINA() function. For example, print() and summary() can be used to display some summary information of model estimation, and AIC() and BIC() can be used to calculate the Akaike (1974) information criterion and the Schwarz (1978) information criterion, respectively. In addition, coef() can be used to extract structural parameters while personparm() can be used to estimate person parameters.
In addition to model estimation, several statistical procedures are available in the GDINA package. First, the Q-matrix can be empirically validated using the Qval() function, which implements de la Torre and Chiu's (2016) ς 2 approach, the ς 2 approach based on empirical cutoffs (Najera et al. 2019), and the stepwise method based on the Wald test (Ma and de la Torre 2020a). A mesa plot (de la Torre and Ma 2016; Ma 2019a) can be created using the plot() function, which is a line graph based on the proportion of variance accounted for by each candidate q-vector (de la Torre and Chiu 2016) to provide a way of visually pinpointing the best q-vector(s) for each item. Second, the modelcomp() function can be used to evaluate whether, for an item which requires at least two attributes, the G-DINA model can be replaced by a reduced model without a significant loss of model-data fit using the Wald test ( (2013), which provide more details about the absolute fit for item pairs, and may be used to identify the sources of misfit. To compare nested models at the test level, an LR test can be conducted using anova(). Additionally, dif() detects differential item functioning using the Wald test (Hou, de la Torre, and Nandakumar 2014) or LR test (Ma et al. 2017).

Illustrations
The data for this illustration were collected from a learning experiment at the University of Tuebingen in Germany in 2010, and previously used by Philipp, Strobl, de la Torre, and Zeileis (2018). Twelve items in elementary probability theory were presented to participants and four attributes were involved: (α 1 ) calculate the probability of the complement of an event, (α 2 ) calculate the probability of two independent events, (α 3 ) calculate the classic probability of an event, and (α 4 ) calculate the probability of the union of two disjoint events. Responses of 504 participants from the first part of the experiment were used. The data and Q-matrix are available from the R package pks (Heller and Wickelmaier 2013). Specifically, after installing the pks package, the data can be loaded by R> data("probability", package = "pks") R> pb <-probability[, sprintf("b1%.2i", 1:12)] Note that the data contain some missing values, which are coded as NA. The Q-matrix is R> Q <-read.table(header = TRUE, text = " + cp id pb un + 0 0 1 0 + 1 0 0 0 + 0 0 0 1 + 0 1 0 0 + 1 0 1 0 + 1 0 1 0 + 0 0 1 1 + 0 0 1 1 + 0 1 1 0 + 1 1 0 0 + 1 1 1 0 + 0 1 1 1") We first fit the G-DINA model to the data by specifying data and Q-matrix as below. A challenge for the EM algorithm is that the solutions could be local maxima. We set nstarts = 200 in argument control to request the function to evaluate the observed log likelihood based on 200 sets of randomly generated initial values. The best set of initial values is used for the following model calibration.

R> coef(GDINA_est, what = "delta")
By specifying what = "lambda", parameters involved in the joint attribute distribution can be extracted. By default, the saturated model is used and therefore, the mixing proportion parameters are printed: Individuals' attribute patterns can be estimated using maximum likelihood estimation (MLE), expected a posteriori (EAP) or maximum a posteriori (MAP) as methods. To obtain attribute estimates, use the function personparm() and specify the object returned from GDINA() as the input. Similar to coef(), personparm() also has an argument called what, where users can specify the method for person attribute estimation. By default, the EAP method is employed. When using MLE or MAP, personparm() will print an additional column consisting of TRUE or FALSE indicating whether the likelihood or posterior have multiple modes or not. If there are multiple modes, the estimate returned is randomly selected. The following code shows the MAP estimates of the first 6 individuals and it can be found that the posterior distributions for all of these 6 individuals are unimodal.

R> map <-personparm(GDINA_est, what = "MAP") R> map[1:6, ]
A1 A2 A3 A4 multimodes 1 1 1 1 1 FALSE 2 1 1 1 1 FALSE 3 1 0 1 1 FALSE 4 1 1 1 1 FALSE 5 1 0 1 1 FALSE 6 1 1 1 1 FALSE To fit other CDMs, we can modify argument model of the GDINA() function. It can be a character vector for each item or a scalar which will be used for all items. For example, to fit the DINA model, we set model = "DINA". We can also modify att.dist to specify the model for the joint attribute distribution. The code below estimates a higher-order DINA model using the Rasch model. Note that the argument higher.order allows further specifications of the higher-order model.
R> HoDINA_est <-GDINA(dat = pb, Q = Q, model = "DINA", + control = list(nstarts = 200), att.dist = "higher.order", + higher.order = list(model = "Rasch")) We can still use summary() to print summary information, coef() to print model parameters involved in item response function and joint attribute distribution, and personparm() to estimate person parmeters (including higher-order ability). For example, to print the structural parameters of the joint attribute distribution, In addition, we can explore whether the G-DINA model can be constrained to some reduced models for each item. To achieve this goal, modelcomp() can be used with the object returned from GDINA() as the input: R> mc <-modelcomp(GDINA.obj = GDINA_est) By default, modelcomp() performs the Wald test for each item that requires two or more attributes. The following code prints the primary results of model comparsions. The column labeled "models" gives the suggested CDM for each item based on the "simpler model + largest p value rule" . Specifically, if the DINA or DINO model was one of the retained models, then the DINA or DINO model with the larger p value was selected as the best model; but if both DINA and DINO were rejected, the reduced model with the largest p value was selected as the best model for this item. Note that when the p values of several reduced models were greater than the nominal level, the DINA and DINO models were preferred over the A-CDM, LLM, and R-RUM models because of their simplicity. The suggested models for the first four items are the G-DINA model because they are all singleattribute items.

R> mc
Item-level model selection: Users can specify method = "LM" or method = "LR" in modelcomp() for score and LR tests. Note that the null hypothesis is that the reduced model can fit the data as well as the G-DINA model, and hence, a non-significant result implies that the reduced model may be used instead of the G-DINA model. Note that modelcomp() is only suitable for nested models: the saturated G-DINA model and the models it subsumes. To compare non-nested models, users can use AIC() and BIC() functions.
Another important function in the GDINA package is Qval() for Q-matrix validation. Based on the object returned from GDINA(), the following code examines whether any element in the Q-matrix is potentially misspecified using the stepwise Wald test (Ma and de la Torre 2020a).
It can be observed that based on the data, a few modifications are suggested to the q-vectors of items 9, 10, 11 and 12. To better understand the results, we can draw mesa plots for these items. For example, the code below draws a mesa plot for item 12 (see Figure 1):

R> plot(Qv, item = 12)
The mesa plot in Figure 1, by default, has the best q-vectors given the number of attributes required on the x-axis. 0 is not a valid q-vector but presented for reference. de la Torre and Ma (2016) noted that the q-vector on the edge of the mesa is potentially the best. In the figure, the one with red dot is the original q-vector. The plot shows that attribute 2 itself (i.e., q j = [0100] ) can explain most variation in the success probabilities, and that attributes 3 and 4 do not contribute much. This implies that attributes 3 and 4 may not be necessary to correctly answer this item. Note that the Q-matrix validation procedures in the GDINA package are based on the data at hand, and whether the suggested modifications should be incorporated should be subject to the judgement of domain experts.

Program comparisons
In this section, we use the R package GDINA to replicate published results based on two sets of real data analyzed using commercial software programs. The first replication is based on 536 students' responses to 20 fraction subtraction items, which has been previously analyzed by many researchers (e.g., de la Torre and Douglas 2004). flexMIRT (Cai 2017) has been used to fit the 1PL higher-order DINA model to the data, and the corresponding syntax and outputs can be found in its manual (Houts and Cai 2016, p. 156-159). The data and Q-matrix can be loaded by: R> data("frac20", package = "GDINA") R> frac_data <-frac20$dat R> frac_Q <-frac20$Q Using the GDINA R package, the 1PL higher-order DINA model can be estimated using the code below: R> est <-GDINA(dat = frac_data, Q = frac_Q, model = "DINA", + att.dist = "higher.order", higher.order = list(model = "1PL", + InterceptRange = c(-5, 5), nquad = 49), + control = list(conv.crit = 1e-6))  Table 2: Higher-order structural parameter estimates from flexMIRT and GDINA.
Note that we specify the 1PL model as the higher-order model by setting model = "1PL" in argument higher.order. We also set the range of intercept parameters at [−5, 5] and the number of quadrature nodes at 49. In addition, we use a more stringent convergence criterion, that is, 1e-6, instead of the default 1e-4. The following code can be used to extract item and higher-order structural parameters: R> gs <-coef(est, what = "gs") R> ho <-coef(est, what = "lambda") The results are given in Tables 1 and 2. The item parameter estimates of flexMIRT in Table 1 Item Templin and Hoffman (2013  were obtained from its manual (Houts and Cai 2016, p. 159), where numbers were reported to only two decimal places. All item parameter estimates in Table 1 from flexMIRT and GDINA are virtually identical. The higher-order structural parameters of flexMIRT in Table 2 were obtained by re-running the flexMIRT code in its manual (see Example 7-7; Houts and Cai 2016, p. 156). The outputs are provided as supplementary material. It can be observed that the GDINA package and flexMIRT also yielded virtually identical estimates of higher-order structural parameter. It is worth noting that the GDINA package completed its calibration in 4.11 seconds, whereas on the same laptop, the flexMIRT calculation took 23.27 seconds.
The second replication is based on the Examination for the Certificate of Proficiency in English (ECPE) data set, which consists of 2922 students' responses to 28 items. The data and Q-matrix can be loaded by: R> data("ecpe", package = "GDINA") R> ecpe_data <-ecpe$dat R> ecpe_Q <-ecpe$Q Templin and Hoffman (2013) fit the loglinear CDM (LCDM; Henson et al. 2009) to the data using Mplus (Muthén and Muthén 2017). Templin and Hoffman (2013) also imposed monotonic constraints to the success probabilities so that mastering of additional attributes will not lead to a lower probability of success. More formally, the monotonic constraint implies that P(α * lj ) ≤ P(α * l j ) if α * lj ≺ α * l j . Since the LCDM is equivalent to the logit link G-DINA model, it can be estimated using the code below by specifying model = "logitGDINA": R> ecpe_est <-GDINA(dat = ecpe_data, Q = ecpe_Q, model = "logitGDINA", + mono.constraint = TRUE, control = list(conv.crit = 1e-6)) In the code above, mono.constraint = TRUE ensures that the monotonic constraints are met during the calibration. Table 3 gives the δ parameter estimates from the GDINA package and Templin and Hoffman (2013) using Mplus. It can be found that the estimates are virtually identical with only a few differences at the third decimal places. We also provided Mplus code as supplemental material, which yielded identical parameter estimates as those in Table 3 from Templin and Hoffman (2013), with only a few exceptions at the third decimal places for items 16 and 23. Note that the GDINA package completed its calibration in 21.42 seconds, whereas Mplus took 32 minutes and 4 seconds.

Summary and discussion
This paper introduces the R package GDINA for conducting a variety of CDM analyses under the G-DINA model framework. There are several features that are potentially useful but were not demonstrated in this paper. First, the G-DINA model has been extended to accommodate polytomous attributes (Chen and de la Torre 2013), polytomous responses (Ma and de la Torre 2016) and multiple-group analysis (Ma et al. 2017). These models can also be estimated using the GDINA() function. Second, the package can also calibrate various models for multiple strategies, including the multiple strategy DINA model (Huo and de la Torre 2014), generalized multiple strategy CDMs (Ma and Guo 2019), and diagnostic tree models (Ma 2019b). Other models and approaches for cognitive diagnosis that can be handled by the package include the multiple-choice DINA model (de la Torre 2009) and iterative latent class analysis (Jiang 2019). Besides, autoGDINA() is a wrapper function to simplify the CDM analyses by conducting the Q-matrix validation, item-level model selection and model calibration sequentially in a one run. In addition, to simulate item responses based on various CDMs under the G-DINA model framework, simGDINA() can be employed. Last but not least, to make the package more accessible to users who are less familar with R, a graphical user interface (GUI), which is developed using shiny (Chang, Cheng, Allaire, Xie, and McPherson 2020) and shinydashboard (Chang and Borges Ribeiro 2018), can be started for various CDM analyses by calling startGDINA(). The GUI session will be launched in the user's default browser. Figure 2 shows the start-up page of the interface, where users could browse and import the item response matrix and Q-matrix in various formats. Figure 3 shows the second tab, where users could choose a CDM from a dropdown list for model calibration or specify different CDMs for different items. Users can also choose the options for validating Q-matrix and performing item-level model selection using the Wald test. After clicking the "CLICK  TO ESTIMATE!" button, a few other tabs will become visible to give various outputs. More details on the use of the GUI for CDM analyses can be found in Ma and de la Torre (2019) and de la Torre and Akbay (2019).
Because of the various functionality, the GDINA package could provide a set of useful tools for researchers and practitioners who are interested in CDMs. The package is being developed actively, and the features that may be incorporated into the package in the future include (1) accommodating complex sampling design, (2) providing other methods for estimating standard errors in the EM algorithm, (3) providing procedures for detecting model identifiability and (4) accommodating a large dataset and Q-matrix.

Computational details
The results in this paper were obtained using R 4.0.0 with the GDINA 2.8.0 package. R itself and all packages used are available from CRAN at https://CRAN.R-project.org/. All analyses were conducted on a laptop with Intel i7-7600U CPU, 8GB RAM and Windows 10 OS.