Cognitively Diagnostic Analysis in R
Introduction
Cognitive diagnosis aims to identify what students have mastered, and what they have yet to.
Cognitive diagnosis models (CDMs) are
latent variable models for cognitive diagnosis
formal models that connect the observed variables (items) with unobserved variables (attributes)
multidimensional in nature to handle multiple latent variables simultaneously
parameterized based on certain “theories” or assumptions
Also called Diagnostic classification models, Cognitive psychometric models, Multiple classification (latent class) models, Latent response models, Restricted latent class models, Structured located latent class models, Structured IRT models
Attributes are
any procedures, skills, or processes that involve in the problem-solving process in cognitive tests
more generally, psychological characteristics (e.g., disorder, personality, attitude) that may affect individuals’ responses to test items
assumed to be discrete (binary most often) latent variables in the CDM context
A toy example
Suppose a math test is created to measure \(K=3\) attributes: addition, subtraction and multiplication.
Each of them is assumed to be binary
1 - mastery
0 - nonmastery
Below are three example items:
3x4-2
6-1
5+3x6
The primary goal of CDMs is to
estimate individuals’ attribute profiles [addition, subtraction, multiplication] or
(equivalently) classify individuals into latent class \(1,2,\dots,2^K=8\) when all profiles are permissible
To ensure the reliability and validity of the classifications, we need to make sure
items have good psychometric properties (item parameter examination, DIF evaluation, etc.)
items measure what they are designed to measure (Q-matrix validation)
make sure CDMs used can fit data adequately (model data fit evaluation)
classifications are reliable (classification reliability)
In what follows, we will first discuss how to use GDINA R package for CDM estimation and then briefly discuss some other analyses when time permits.
CDMs for binary data and binary attributes
There are a large number of CDMs in the literature (Davier and Lee 2019). We are going to briefly review a few of them here.
Suppose there are \(K\) binary attributes, we have \(2^K\) attribute profiles, each labelling a unique latent class. The attribute profile for latent class \(c\) is denoted as \(\boldsymbol{\alpha}_c=[\alpha_{c1},\ldots,\alpha_{cK}]^\top\), where \(\alpha_{ck}=1\) if attribute \(k\) is mastered and \(\alpha_{ck}=0\) if not.
For the toy example above, there are 8 latent classes:
| Latent Class | Addition | Subtraction | Multiplication |
|---|---|---|---|
| 1 | 0 | 0 | 0 |
| 2 | 1 | 0 | 0 |
| 3 | 0 | 1 | 0 |
| 4 | 0 | 0 | 1 |
| 5 | 1 | 1 | 0 |
| 6 | 1 | 0 | 1 |
| 7 | 0 | 1 | 1 |
| 8 | 1 | 1 | 1 |
However, if an item does not measure all attributes, some latent classes have the same success probabilities on this item. For example, Item 1 does not measure addition, so mastering addition is not supposed to affect the success probabilities on this item.
Some latent classes can be collapsed and we call them latent groups.
For item 1 again, we have the following four latent groups. With De La Torre (2011) ’s notation, those (sub-)profiles are denoted by \(\boldsymbol{\alpha}_{lj}^*\) for item \(j\) and latent group \(l\).
| Latent Group | Subtraction | Multiplication |
|---|---|---|
| 1 | 0 | 0 |
| 2 | 1 | 0 |
| 3 | 0 | 1 |
| 4 | 1 | 1 |
The DINA model
For the DINA model (Haertel 1989), the item response function can be written by \[ P(X_j=1|\boldsymbol{\alpha}_{lj}^*) = \delta_{j0} + \delta_{j1} \prod_{k=1}^{K_j} \alpha_{lk} \]
The G-DINA model
The G-DINA model is a generalized DINA model (De La Torre 2011), where different latent groups can have distinct success probabilities. \[ P(X_j=1|\boldsymbol{\alpha}_{lj}^*) =\delta_{j0}+\sum_{k=1}^{K_j^*}\delta_{jk}\alpha_{lk}+\sum_{k'=k+1}^{K_j^*}\sum_{k=1}^{K_j^*-1}\delta_{jkk'}\alpha_{lk}\alpha_{lk'}+\cdots+\delta_{j12{\cdots}K_j^*}\prod_{k=1}^{K_j^*}\alpha_{lk}, \]
CDM Estimation in R
Installation
Please download and install the latest R software (version 4.4.1) from https://cran.r-project.org/.
In addition, please download and install the latest RStudio, which is a free graphic user interface for R and will be used in this training session. It can be downloaded from https://www.rstudio.com
After you install R and RStudio, open RStudio and install the GDINA R package:
In the Packages panel, you should be able to see the GDINA package version 2.9.4:
Estimation of the G-DINA Model
data1 consists of responses of 837 individuals to 15 items, generated based on a proportional reasoning assessment. Q1 is the associated Q-matrix, involving 3 attributes.
In this workshop, I will use data1 for illustration.
# data
data1 <- read.table(file = "data1.dat", header = TRUE)
head(data1) Item1 Item2 Item3 Item4 Item5 Item6 Item7 Item8 Item9 Item10 Item11 Item12
1 1 1 1 1 1 0 1 1 0 1 1 1
2 1 1 1 1 0 0 0 0 1 0 1 1
3 1 0 1 0 1 1 1 1 1 0 1 0
4 1 0 1 0 1 0 0 0 1 1 1 0
5 1 0 0 0 0 0 0 0 0 0 1 0
6 1 0 1 1 1 1 1 1 0 1 1 0
Item13 Item14 Item15
1 0 1 0
2 0 0 0
3 0 0 0
4 0 1 0
5 0 1 1
6 1 0 0# Q-matrix
Q1 <- read.table(file = "Q1.txt")
Q1 V1 V2 V3
1 1 0 0
2 0 0 1
3 0 1 0
4 0 0 1
5 1 1 0
6 1 1 0
7 1 1 0
8 0 1 1
9 0 1 1
10 0 1 1
11 0 1 1
12 1 0 1
13 1 0 1
14 1 1 1
15 1 1 1The GDINA function is the main function of the package. It takes various arguments as in:
R Code
GDINA(dat, Q, model = "GDINA", ...)
To estimate G-DINA model, call GDINA function and specify the data and Q-matrix as the first two arguments.
# Load the package
library(GDINA)GDINA R Package (version 2.9.4; 2023-07-01)
For tutorials, see https://wenchao-ma.github.io/GDINA# Fit the data using G-DINA model
fit1 <- GDINA(dat = data1, Q = Q1, model = "GDINA")Returned fit1 contains (almost) all estimation results. It is called an object of class GDINA, and we can apply various methods to it.
Methods for printing summary information
To print some general model estimation information, type fit1 in Rstudio console:
fit1Call:
GDINA(dat = data1, Q = Q1, model = "GDINA")
GDINA version 2.9.4 (2023-07-01)
===============================================
Data
-----------------------------------------------
# of individuals groups items
837 1 15
===============================================
Model
-----------------------------------------------
Fitted model(s) = GDINA
Attribute structure = saturated
Attribute level = Dichotomous
===============================================
Estimation
-----------------------------------------------
Number of iterations = 288
For the final iteration:
Max abs change in item success prob. = 0.0001
Max abs change in mixing proportions = 0.0000
Change in -2 log-likelihood = 0.0001
Converged? = TRUE
Time used = 0.3014 secs As shown above, some information about the data, fitted model and calibration was printed. summary is another method that can be applied to objects of class GDINA. Specify fit1 as the input:
summary(fit1)
Test Fit Statistics
Loglik = -7431.45
AIC = 14996.91 | penalty [2 * p] = 134.00
BIC = 15313.81 | penalty [log(n) * p] = 450.90
CAIC = 15380.81 | penalty [(log(n) + 1) * p] = 517.90
SABIC = 15101.04 | penalty [log((n + 2)/24) * p] = 238.13
No. of parameters (p) = 67
No. of estimated item parameters = 60
No. of fixed item parameters = 0
No. of distribution parameters = 7
Attribute Prevalence
Level0 Level1
A1 0.3688 0.6312
A2 0.4601 0.5399
A3 0.2875 0.7125Methods for item parameters
R code
coef(object, what = c(“catprob”, “delta”, “gs”, “lambda”, …), withSE,…)
To extract estimated item parameters \(\hat{\boldsymbol{\delta}}\)s, we can use coef function with what = "delta":
coef(object = fit1, what = "delta")$`Item 1`
d0 d1
0.9291 0.0517
$`Item 2`
d0 d1
0.2614 0.5588
$`Item 3`
d0 d1
0.5381 0.2868
$`Item 4`
d0 d1
0.4810 0.3897
$`Item 5`
d0 d1 d2 d12
0.2881 0.4402 0.2551 -0.0912
$`Item 6`
d0 d1 d2 d12
0.1571 0.2843 -0.1107 0.4262
$`Item 7`
d0 d1 d2 d12
0.0001 0.5171 0.5315 -0.2725
$`Item 8`
d0 d1 d2 d12
0.0395 0.6293 0.3077 -0.1728
$`Item 9`
d0 d1 d2 d12
0.0937 0.3138 0.2515 0.0233
$`Item 10`
d0 d1 d2 d12
0.4337 0.2776 -0.0742 0.1376
$`Item 11`
d0 d1 d2 d12
0.4877 0.1327 0.2384 0.0473
$`Item 12`
d0 d1 d2 d12
0.1883 0.1551 0.1680 0.1733
$`Item 13`
d0 d1 d2 d12
0.1980 -0.0068 0.1097 0.2900
$`Item 14`
d0 d1 d2 d3 d12 d13 d23 d123
0.0764 0.9235 -0.0523 0.1596 -0.3355 -0.7489 -0.0965 0.7242
$`Item 15`
d0 d1 d2 d3 d12 d13 d23 d123
0.0989 0.6058 -0.0917 0.2013 -0.2809 -0.5753 0.0388 0.6482 We are often also interested in item success probabilities. To obtain \(\hat{P}(X_j=1|\boldsymbol{\alpha}_{lj}^*)\), use coef but specify what = "catprob". If you would like to print standard errors of item parameters, simply set withSE=TRUE.
coef(object = fit1, what = "catprob", withSE = TRUE)$`Item 1`
P(0) P(1)
Est. 0.9291 0.9808
S.E. 0.0172 0.0073
$`Item 2`
P(0) P(1)
Est. 0.2613 0.8201
S.E. 0.0486 0.0224
$`Item 3`
P(0) P(1)
Est. 0.5381 0.8249
S.E. 0.0309 0.0228
$`Item 4`
P(0) P(1)
Est. 0.4810 0.8707
S.E. 0.0456 0.0183
$`Item 5`
P(00) P(10) P(01) P(11)
Est. 0.2881 0.7283 0.5432 0.8921
S.E. 0.0448 0.0567 0.0726 0.0233
$`Item 6`
P(00) P(10) P(01) P(11)
Est. 0.1571 0.4414 0.0464 0.7568
S.E. 0.0359 0.0598 0.0768 0.0332
$`Item 7`
P(00) P(10) P(01) P(11)
Est. 0.0001 0.5172 0.5316 0.7762
S.E. 0.0369 0.0658 0.0774 0.0291
$`Item 8`
P(00) P(10) P(01) P(11)
Est. 0.0395 0.6688 0.3473 0.8037
S.E. 0.0746 0.0539 0.0375 0.0347
$`Item 9`
P(00) P(10) P(01) P(11)
Est. 0.0937 0.4075 0.3453 0.6824
S.E. 0.0676 0.0536 0.0371 0.0375
$`Item 10`
P(00) P(10) P(01) P(11)
Est. 0.4337 0.7113 0.3595 0.7748
S.E. 0.0844 0.0502 0.0383 0.0355
$`Item 11`
P(00) P(10) P(01) P(11)
Est. 0.4877 0.6204 0.7261 0.9061
S.E. 0.0855 0.0509 0.0343 0.0258
$`Item 12`
P(00) P(10) P(01) P(11)
Est. 0.1883 0.3434 0.3564 0.6847
S.E. 0.0676 0.0529 0.0445 0.0319
$`Item 13`
P(00) P(10) P(01) P(11)
Est. 0.1980 0.1912 0.3077 0.5909
S.E. 0.0694 0.0474 0.0433 0.0332
$`Item 14`
P(000) P(100) P(010) P(001) P(110) P(101) P(011) P(111)
Est. 0.0764 0.9999 0.0241 0.2360 0.6121 0.4106 0.0872 0.6505
S.E. 0.0782 0.2331 0.1439 0.0509 0.0637 0.0666 0.0804 0.0444
$`Item 15`
P(000) P(100) P(010) P(001) P(110) P(101) P(011) P(111)
Est. 0.0989 0.7047 0.0072 0.3002 0.3321 0.3307 0.2473 0.6451
S.E. 0.0856 0.1520 0.1438 0.0548 0.0602 0.0640 0.0939 0.0452If you want to print the success probabilities for all latent classes, \(\hat{P}(X_j=1|\boldsymbol{\alpha}_{c})\), you can specify what = "LCprob".
coef(object = fit1, what = "LCprob") 000 100 010 001 110 101 011 111
Item 1 0.9291 0.9808 0.9291 0.9291 0.9808 0.9808 0.9291 0.9808
Item 2 0.2613 0.2613 0.2613 0.8201 0.2613 0.8201 0.8201 0.8201
Item 3 0.5381 0.5381 0.8249 0.5381 0.8249 0.5381 0.8249 0.8249
Item 4 0.4810 0.4810 0.4810 0.8707 0.4810 0.8707 0.8707 0.8707
Item 5 0.2881 0.7283 0.5432 0.2881 0.8921 0.7283 0.5432 0.8921
Item 6 0.1571 0.4414 0.0464 0.1571 0.7568 0.4414 0.0464 0.7568
Item 7 0.0001 0.5172 0.5316 0.0001 0.7762 0.5172 0.5316 0.7762
Item 8 0.0395 0.0395 0.6688 0.3473 0.6688 0.3473 0.8037 0.8037
Item 9 0.0937 0.0937 0.4075 0.3453 0.4075 0.3453 0.6824 0.6824
Item 10 0.4337 0.4337 0.7113 0.3595 0.7113 0.3595 0.7748 0.7748
Item 11 0.4877 0.4877 0.6204 0.7261 0.6204 0.7261 0.9061 0.9061
Item 12 0.1883 0.3434 0.1883 0.3564 0.3434 0.6847 0.3564 0.6847
Item 13 0.1980 0.1912 0.1980 0.3077 0.1912 0.5909 0.3077 0.5909
Item 14 0.0764 0.9999 0.0241 0.2360 0.6121 0.4106 0.0872 0.6505
Item 15 0.0989 0.7047 0.0072 0.3002 0.3321 0.3307 0.2473 0.6451If you just want to print the guessing and slip parameters, or \(\hat{P}(X_j=1|\boldsymbol{\alpha}_{c}=\boldsymbol{0})\) and \(\hat{P}(X_j=0|\boldsymbol{\alpha}_{c}=\boldsymbol{1})\) respectively, you can specify what = "gs".
coef(object = fit1, what = "gs") guessing slip
Item 1 0.9291 0.0192
Item 2 0.2613 0.1799
Item 3 0.5381 0.1751
Item 4 0.4810 0.1293
Item 5 0.2881 0.1079
Item 6 0.1571 0.2432
Item 7 0.0001 0.2238
Item 8 0.0395 0.1963
Item 9 0.0937 0.3176
Item 10 0.4337 0.2252
Item 11 0.4877 0.0939
Item 12 0.1883 0.3153
Item 13 0.1980 0.4091
Item 14 0.0764 0.3495
Item 15 0.0989 0.3549Methods for person parameters
R code
personparm(fit1, what = c(“EAP”, “MAP”, “MLE”, “mp”) )
By default, personparm function estimates attribute status using EAP method.
#EAP estimates
head(personparm(fit1)) A1 A2 A3
[1,] 1 1 1
[2,] 0 0 1
[3,] 1 1 0
[4,] 1 1 0
[5,] 1 0 0
[6,] 1 1 0You can however, use MAP or MLE. One additional column is printed to indicate whether there are multiple modes in posterior and likelihood.
#MAP estimates
head(personparm(fit1, what = "MAP")) A1 A2 A3 multimodes
1 1 1 1 FALSE
2 0 0 1 FALSE
3 1 1 0 FALSE
4 1 1 0 FALSE
5 1 0 0 FALSE
6 1 1 0 FALSE#MLE estimates
head(personparm(fit1, what = "MLE")) A1 A2 A3 multimodes
1 1 1 1 FALSE
2 0 0 1 FALSE
3 1 1 0 FALSE
4 1 1 0 FALSE
5 1 0 0 FALSE
6 1 1 0 FALSEWhen specifying what = "mp", you will get the mastery probabilities of each individual to each attribute.
#probability of mastering each attribute
head(personparm(fit1, what = "mp")) A1 A2 A3
[1,] 0.9174 0.8540 0.8797
[2,] 0.0889 0.0709 0.9849
[3,] 0.9714 0.9942 0.0348
[4,] 0.8143 0.7877 0.1100
[5,] 0.7215 0.0048 0.1471
[6,] 0.9688 0.9516 0.3836Methods for plots and extracting other information
R code
plot(x, what = c(“IRF”, “mp”), item, person,…)
extract(object, what,…)
To plot the probabilities of success for each item \(\hat{P}(X_{j}=1|\boldsymbol{\alpha}_{lj}^*)\), specify what = "IRF":
plot(fit1, what = "IRF", item = c(9,15), withSE=TRUE)
To plot the mastery probability, specify what = "mp":
plot(fit1, what = "mp", person = c(1,5))
To extract other elements, we should use extract function. Its first input should be an object returned from GDINA function (or other functions), and its second input specifies what to extract. See its help page ?extract for a complete list of elements that can be extracted. For example, the posterior probabilities can be printed:
extract(fit1, what = "posterior.prob") 000 100 010 001 110 101 011
[1,] 0.06151017 0.03144658 0.03741109 0.186019 0.1570853 0.1811396 0.08382109
111
[1,] 0.2615672Some other considerations in model estimation
Monotonic constraints
higher-order structure for attributes
attributes hierarchical structure
A more complete CDM analysis
Shi et al. (2021) created a CDM analysis flowchart:
Q-matrix validation
GDINA package has several implementations of Q-matrix validation. Below is an example:
Qval(fit1, method = "wald")
Q-matrix validation based on Stepwise Wald test
Suggested Q-matrix:
A1 A2 A3
1 1 0 0
2 0 0 1
3 0 1 0
4 0 0 1
5 1 1 0
6 1 1 0
7 1 1 0
8 0 1 1
9 0 1 1
10 0 1 0*
11 0 1 1
12 1 0 1
13 1 0 1
14 1 0* 0*
15 1 1 1
Note: * denotes a modified element.Model fit evaluation
To compare the absolute fit of a model to the data, one can use modelfit and itemfit:
modelfit(fit1)Test-level Model Fit Evaluation
Relative fit statistics:
-2 log likelihood = 14862.91 ( number of parameters = 67 )
AIC = 14996.91 BIC = 15313.81
CAIC = 15380.81 SABIC = 15101.04
Absolute fit statistics:
M2 = 60.8605 df = 53 p = 0.214
RMSEA2 = 0.0133 with 90 % CI: [ 0 , 0.0266 ]
SRMSR = 0.0269itemfit allows us to evaluate the fit between item pairs:
ift <- itemfit(fit1)
plot(ift)
Exercise
Please analyze the ECPE data, which consists of responses of 2922 examinees to 28 items involving 3 attributes. Attribute 1 is morphosyntactic rules, Attribute 2 is cohesive rules and Attribute 3 is lexical rules.
data <- GDINA::ecpe$dat
Q <- GDINA::ecpe$QConduct Q-matrix validation and see if any modifications are recommended.
With the original Q-matrix, fit both G-DINA model and DINA model to the data.
compare two models via
anovafunction.evaluate the absolute fit of the preferred model.
Find the proportion of students who master each attribute based on preferred model.
Find the proportion of students with each attribute profile.
Appendix: Additional Information
What can the GDINA package do?
- Estimating G-DINA model and a variety of widely-used models subsumed by the G-DINA model, including the DINA model, DINO model, additive-CDM (A-CDM), linear logistic model (LLM), reduced reparametrized unified model (RRUM), multiple-strategy DINA model for dichotomous responses
- Estimating models within the G-DINA model framework using user-specified design matrix and link functions
- Estimating Bugs-DINA, DINO and G-DINA models for dichotomous responses
- Estimating sequential G-DINA model for ordinal and nominal responses
- Estimating polytomous G-DINA model for ordinal attributes
- Estimating generalized multiple-strategy CDM and diagnostic tree model for multiple strategies
- Estimating an extended multiple-choice DINA model
- Modeling independent, saturated, higher-order, loglinear smoothed, and structured joint attribute distribution
- Accommodating multiple-group analysis
- Imposing monotonic constrained success probabilities
- Validating Q-matrix under the general model framework using various approaches
- Evaluating absolute and relative item and model fit
- Comparing models at the test and item levels
- Detecting differential item functioning using Wald and likelihood ratio test
- Evaluating classification accuracy
- Simulating data based on all aforementioned CDMs
- Providing graphical user interface for users less familiar with R
More Learning Resources
Package website with many additional examples and tutorials: https://wenchao-ma.github.io/GDINA
General introductions:
“Official” introduction to the package published in Journal of Statistical Software: https://doi.org/10.18637/jss.v093.i14
A book chapter on an introduction to the GDINA package: https://doi.org/10.1007/978-3-030-05584-4_29
On the use of graphical user interface:
An article on how to conduct various CDM analyses using the graphical user interface: https://doi.org/10.14689/ejer.2019.80.9
An NCME digital module on a gentle introduction to the G-DINA model framework and the use of graphical user interface for CDM analyses: https://ncme.elevate.commpartners.com/
References
Davier, Matthias von, and Young-Sun Lee. 2019. “Handbook of Diagnostic Classification Models.” Cham: Springer International Publishing.
De La Torre, Jimmy. 2011. “The Generalized DINA Model Framework.” Psychometrika 76: 179–99.
Haertel, Edward H. 1989. “Using Restricted Latent Class Models to Map the Skill Structure of Achievement Items.” Journal of Educational Measurement 26 (4): 301–21.
Shi, Qingzhou, Wenchao Ma, Alexander Robitzsch, Miguel A Sorrel, and Kaiwen Man. 2021. “Cognitively Diagnostic Analysis Using the g-DINA Model in r.” Psych 3 (4): 812–35.