Latent classes and cognitive diagnosis

Module: Beyond the standard models
Tranche: extension (beyond a first course)
Before this: What is an item? Explanatory item response models
Builds toward: none
Picks up threads: The sum score is sufficient under the Rasch model (from The Rasch model); Items are a sample from a population of possible items, describable by their features (from What is an item? Explanatory item response models); Ordered latent classes can look like a continuum cut into pieces (from Constructs and construct maps)
IRW tables: frac20, cdm_ecpe
From EDUC 252: slides c7; problem sets PS7#2; code c7/3_cdm.R, ps7/cdm.R
Status: draft
Authorship: Written largely by Claude (Anthropic), from Ben Domingue’s EDUC 252 materials and under Ben’s direction. Ben reviews each lesson before it is marked done.
Comments: open an issue on GitHub about this lesson.

What this is for

What is an item? explained item difficulty with a matrix of item features and kept every respondent on one continuum. This lesson keeps the matrix and drops the continuum. Suppose respondents differ not in how much of something they have but in which skills they have: a student who can find a common denominator but can’t borrow from a whole number needs different teaching from one with the opposite gap, even if both get half the items right. Latent class analysis puts respondents in a few unobserved groups; cognitive diagnosis builds the groups from profiles of skills, with a matrix saying which skills each item needs. Constructs and construct maps asked whether a construct is a continuum or a set of classes, and summary statistics couldn’t tell. Here we fit both and see which the data prefer.

Goals

By the end of this lesson you should be able to:

  • write down the latent class model, fit it with poLCA, and choose the number of classes;
  • say why ordered latent classes behave like a continuum cut into pieces, and what would be evidence of something else;
  • read a Q-matrix and say which items a skill profile should answer correctly under DINA;
  • fit DINA and G-DINA with GDINA, and read guessing, slipping and the share of respondents with each skill;
  • compare a cognitive diagnosis model with the Rasch model and the 2PL, and use empirical Q-matrix validation with care.

Core ideas

Latent classes

Recall from Constructs and construct maps: latent classes are unobserved groups inside each of which people are alike. Ordered classes close together gave nearly the same sum scores and item-rest correlations as a continuum, and telling them apart was left to models fit to the responses.

A latent class model says each respondent belongs to one of \(C\) classes, class \(c\) holding a share \(\pi_c\) of the population, and that within a class each item has its own probability of success, \(p_{ic}\). Given the class, responses are independent: local independence, with a class in place of \(\theta\). The probability of respondent \(j\)’s responses is a mixture over the classes, \[\Pr(\mathbf{x}_j) = \sum_{c=1}^{C} \pi_c \prod_{i} p_{ic}^{\,x_{ij}}\,(1 - p_{ic})^{1 - x_{ij}}.\] The model comes from Lazarsfeld’s latent structure analysis (Lazarsfeld & Henry, 1968); Goodman (1974) took up its estimation and when its parameters are identified. With \(I\) items it has \(C - 1\) class shares and \(C \times I\) item probabilities. It is fit by EM: from a guess at the item probabilities, compute each respondent’s posterior probability of each class, re-estimate the item probabilities with responses weighted by those posteriors, and repeat. We use the poLCA package (Linzer & Lewis, 2011).

Read the formula from the inside out. Inside class \(c\), each item contributes \(p_{ic}\) if it was answered correctly and \(1 - p_{ic}\) if not, and local independence lets us multiply those factors into the probability of the whole pattern if the respondent were in class \(c\). We don’t know the class, so each class’s product is weighted by its share \(\pi_c\) and the weighted products are added. Each term’s part of the sum is the posterior probability that the respondent is in that class, which is how EM assigns respondents.

Try it: building the probability of a response pattern

Three classes answer four items. Low and High are ordered; Uneven finds items 1 and 2 easy and items 3 and 4 hard. Choose a response pattern and the classes’ relative sizes.

With the default pattern, items 1 and 2 right and 3 and 4 wrong, Uneven supplies most of the probability (a posterior of 0.63), though it is the smallest class. Now answer items 3 and 4 instead: the score is still 2, but the pattern is about thirty times less likely and Low takes 0.70 of the posterior. A latent class model reads which items were answered, not only how many. Setting a class’s size to 0 removes its term from the sum.

If you’ve done Estimating item parameters: the two biased coins that introduced EM were a two-class latent class model, with the ten tosses in a set as its items.

Another class always raises the log likelihood, so the number of classes is chosen by a criterion that charges for parameters, such as BIC, and by whether the classes mean anything. A latent class model makes no assumption of a continuum, which is what a theory of types, such as Moffitt’s account of antisocial behaviour, calls for; it can also describe how responses hang together without one.

If you’ve done Model fit and out-of-sample prediction: AIC and BIC, which charge for each parameter, also compare models that aren’t nested, such as a latent class model and a 2PL.

Classes or a continuum?

A latent class model has no \(\theta\). All it has is a set of classes, each a share of respondents and a row of success rates, one per item, like the rows of the table in the last widget. So the question “types or a continuum?” has to be asked of those rows. What pattern in them would a continuum leave behind?

Draw each class’s row as a line across the items, easiest to hardest: a profile. If the classes are really types, respondents who differ in which items they find hard, the profiles cross: one class beats another on some items and loses on others. If the classes are levels of one thing, each class does better than the class below it on every item, and the profiles stack without crossing.

Try it: class profiles, stacked or crossing

The widget makes classes that are levels of one thing. The top panel shows where they come from, which the latent class model never sees: a normal \(\theta\) split into equal-sized groups. The bottom panel shows what the model does see, each class’s profile on ten Rasch items, in the same colour as its group. The toggle adds a class that is a type: it finds the odd-numbered items easy and the even-numbered ones hard.

However many levels you choose, the profiles stack: 0 of 10 items are out of order. Add the type and its profile zigzags across the rest, putting all ten out of order. So a latent class model that finds stacked profiles has found what a continuum would give it, a continuum cut into a few levels, not evidence of types. The converse holds too: a Rasch model whose \(\theta\) takes only a few values, ordered classes by another name, gives nearly the same item estimates as the Rasch model’s conditional likelihood once it has at least \((I + 1)/2\) of them (Lindsay, Clogg & Grego, 1991). Evidence for classes as something other than a continuum has to come from crossings, classes that differ in which items they find hard.

Cognitive diagnosis: classes built from attributes

Recall from What is an item?: the LLTM wrote each item’s easiness as \(\beta_i = \sum_k q_{ik}\eta_k\), where \(q_{ik}\) is 1 if item \(i\) has feature \(k\). The matrix of \(q_{ik}\) described the items, and respondents stayed on one continuum.

Cognitive diagnosis uses such a matrix differently. Each respondent has a profile \(\boldsymbol\alpha_j = (\alpha_{j1}, \dots, \alpha_{jK})\) of binary attributes, skills they have (\(\alpha_{jk} = 1\)) or lack. \(K\) attributes give \(2^K\) profiles, and each profile is a latent class. The Q-matrix says which attributes each item requires: \(q_{ik} = 1\) if item \(i\) needs attribute \(k\). In the LLTM the matrix explained how hard an item is; here it decides which classes an item can tell apart. A cognitive diagnosis model (CDM) is therefore a restricted latent class model (Haertel, 1989): the classes are fixed in advance by the attributes, and the Q-matrix ties their item probabilities together (Rupp, Templin & Henson, 2010).

The classic example is fraction subtraction, from Kikumi Tatsuoka’s work on diagnosing students’ errors (Tatsuoka, 1983; the data as analysed by Tatsuoka, 2002). Twenty items are coded on eight attributes, giving \(2^8 = 256\) profiles. The attribute labels are those of the CDM package documentation, which takes the Q-matrix from DeCarlo (2011) and de la Torre and Douglas (2004).

Try it: the Q-matrix at work

Choose a student’s attributes. Under DINA (next section), an item is answered correctly with probability \(1 - s\) if the student has every attribute it needs, and \(g\) otherwise.

The default student lacks only A5, borrowing from the whole number part; eight items need it, so 12 of 20 are within reach. Now tick A5 and untick A7, subtracting numerators: 19 items need A7, and one is left. Profiles that lack different attributes can also reach similar expected scores through different items, a difference a sum score can’t see.

A profile also describes a respondent, so it can be handed to a simulator. Zhou, Liu, Tang and Wilson (2026), a preprint, took 15 of these fraction items coded on five attributes, wrote each attribute profile as a short description of a student, and asked language models to answer as that student. Answering with no profile, the models were more accurate and more alike than the 536 students. Given profiles, their scores spread more like the students’, and the 1PL difficulties from their answers came closer to the students’ difficulties: in the strongest case, the rank correlation rose from 0.24 to about 0.9. AI and psychometrics takes up language models as simulated respondents more generally.

DINA: all or nothing

The DINA model (“deterministic inputs, noisy and gate”; Junker & Sijtsma, 2001) asks one question per item: does the respondent have every attribute it needs? Write \[\eta_{ij} = \prod_{k} \alpha_{jk}^{\,q_{ik}},\] which is 1 if they do and 0 if they lack any. Then \[\Pr(x_{ij} = 1 \mid \boldsymbol\alpha_j) = (1 - s_i)^{\eta_{ij}}\, g_i^{\,1 - \eta_{ij}}.\] The slip \(s_i\) is the chance that a respondent with everything the item needs still gets it wrong; the guess \(g_i\) is the chance that one without gets it right. Each item has two parameters, because the model splits the \(2^K\) classes into two groups per item. Having some of the required attributes counts for nothing.

If you’ve done From the 1PL to the 4PL: its section on ceilings and slips pointed here. DINA’s \(1 - s\) and \(g\) are an upper and a lower asymptote, like the 4PL’s \(u\) and \(c\), with nothing in between.

Other models relax the rule. In DINO, any one required attribute is enough (Templin & Henson, 2006), which suits a symptom that several conditions can produce. G-DINA (de la Torre, 2011) gives each combination of an item’s \(K_i\) required attributes its own probability, \(2^{K_i}\) of them, and contains DINA and DINO; the log-linear CDM (Henson, Templin & Willse, 2009) is an equivalent family written as a logistic regression on the attributes and their interactions. DINA assumes an item is solved by a chain of steps that must all succeed, as fraction subtraction may be, and spends two parameters per item; G-DINA lets the data say how partial mastery counts, at the price of \(2^{K_i}\) parameters per item, which a small sample may not support. How well each predicts held-out responses says which trade pays in a given dataset (problem 3).

Tying the attributes to a continuum: the higher-order model

DINA says how an item depends on a profile, but not how profiles are spread among respondents. Left free, that spread is a share \(\pi\) for every one of the \(2^K\) profiles, with nothing linking them. Do respondents who have one attribute tend to have the others? The higher-order model (de la Torre & Douglas, 2004) answers yes, through one continuous trait. Each respondent has a \(\theta_j\), normally distributed, and given \(\theta_j\) the attributes are independent, with \[\Pr(\alpha_{jk} = 1 \mid \theta_j) = \frac{\exp(\lambda(\theta_j - \delta_k))}{1 + \exp(\lambda(\theta_j - \delta_k))}.\] Each attribute has a threshold \(\delta_k\), the point on \(\theta\) where it is as likely mastered as not, and \(\lambda\) says how sharply mastery rises around it. This is a Rasch-type model with the attributes in place of the items, except that the attributes are never observed. The share of a profile is the product of these probabilities, averaged over \(\theta\), so the \(2^K - 1\) free shares become \(K\) thresholds (and \(\lambda\)). With \(\lambda\) fixed at 1, as GDINA does by default, that is one parameter per attribute.

What does an item look like against \(\theta\) now? Below the attributes’ thresholds a respondent probably lacks something the item needs and sits near the guess; above them they probably have it all and sit near \(1 - s\).

Try it: two plateaus, and a curve between them

Take one item and a respondent with a given \(\theta\). With thresholds \(\delta_k\) spread from −1 to 1, the widget computes the chance that they have every attribute the item needs, \(\prod_k \Pr(\alpha_{jk} = 1 \mid \theta_j)\) over those attributes, and mixes the plateaus in that proportion: \(\Pr(x_{ij} = 1 \mid \theta_j) = g + (1 - s - g)\prod_k \Pr(\alpha_{jk} = 1 \mid \theta_j)\).

The dashed lines are the plateaus, \(g\) and \(1 - s\), where respondents lacking an attribute and those with all of them sit. Given only \(\theta\) the attributes are uncertain, so the curve runs between. With one attribute it is a 4PL curve (see the Go deeper callout). Require more and the curve moves right and loses its symmetry, since every threshold must be passed. Raise \(\lambda\) and the curve steepens toward a step.

How many classes, which Q-matrix, and why not IRT?

Eight attributes give 256 profiles, so an unrestricted distribution over them has 255 parameters, to be estimated from 536 students in the fraction data. Two restrictions shrink that. The higher-order model above needs 8, one per attribute. An attribute hierarchy allows only profiles consistent with an order of learning, ruling out an advanced skill without its prerequisite (Templin & Bradshaw, 2014). Both move a CDM toward a continuum.

Try it: the cost of a profile

At eight attributes and 536 respondents, an even spread leaves about two students per profile. Each attribute doubles the unrestricted count and adds one to the higher-order count.

The Q-matrix usually comes from content experts, and may be wrong. Empirical Q-matrix validation looks for better entries. The PVAF method (de la Torre & Chiu, 2016) asks, item by item, for the smallest set of attributes whose G-DINA success probabilities carry at least a set share (by default 95%) of the variation that the item’s full set of attributes would give. It can point to items whose coding deserves a second look; a suggestion is a hypothesis about the item, to be checked against what the item asks and by refitting. Is a CDM better than IRT? The models aren’t nested, so we compare them by AIC and BIC, or by predicting held-out responses (problem 3). Ma, Liu, Kanopka, Ma and Domingue (2026) did the latter, with the IMV, on 11 datasets previously analysed with CDMs, including the ECPE and the fraction data used below. When each respondent was assigned their most probable profile, the 2PL predicted held-out responses better than DINA and G-DINA in every dataset. When predictions used each respondent’s probability of mastering each attribute, a CDM did better in 6 of the 11, by small margins (an IMV under 0.02). G-DINA in particular fit the data it was estimated on much better than it predicted new responses. Why might a CDM lose on data it was built for? The IMV site has an interactive version of the simulation side: data come from DINA with correlated skills, and as the correlation rises the profiles collapse toward one dimension and the 2PL closes the gap.

Assumptions: the DINA model for item \(i\), with guess \(g\) and slip \(s\); a higher-order model in which, given a continuous \(\theta\), attribute \(k\) is mastered with probability \(\text{logistic}(\lambda(\theta - \delta_k))\), independently across attributes.

Given \(\theta\), a respondent answers correctly either by having every required attribute and not slipping, or by lacking one and guessing: \[\Pr(x_i = 1 \mid \theta) = (1 - s)\,P_{\text{all}}(\theta) + g\,\big(1 - P_{\text{all}}(\theta)\big) = g + (1 - s - g)\,P_{\text{all}}(\theta),\] where \(P_{\text{all}}(\theta) = \prod_{k:\,q_{ik}=1} \text{logistic}(\lambda(\theta - \delta_k))\) is the probability of having all of them.

With one required attribute, \(P_{\text{all}}(\theta) = \text{logistic}(\lambda(\theta - \delta))\), and \[\Pr(x_i = 1 \mid \theta) = g + \big((1 - s) - g\big)\,\frac{\exp(\lambda(\theta - \delta))}{1 + \exp(\lambda(\theta - \delta))},\] which is the 4PL with \(c = g\), \(u = 1 - s\), slope \(a = \lambda\) and difficulty \(b = \delta\). The curve’s location and slope belong to the attribute: every item that needs only attribute \(k\) has the same slope and difficulty and differs only in its asymptotes.

With \(m \ge 2\) attributes, \(P_{\text{all}}\) is a product of logistic curves, which isn’t a logistic curve. In log-odds its lower tail has slope \(m\lambda\) and its upper tail slope \(\lambda\), so it leaves the floor abruptly and approaches the top slowly; and at the largest threshold every other factor is above 0.5 and that one equals 0.5, so its midpoint sits to the right of every \(\delta_k\). As \(\lambda\) grows, each factor becomes a step at \(\delta_k\) and the product a step at the largest threshold. So the higher-order DINA is a unidimensional IRT model whose item curves are built from the Q-matrix. DeCarlo (2011) fits a logistic version of DINA with software for latent class analysis, along the same lines.

Simulate

The code simulates 1,000 respondents answering ten items on three attributes and fits DINA with GDINA (Ma & de la Torre, 2020), then G-DINA, then DINA with one wrong Q-matrix entry. The first run takes about half a minute, most of it installing GDINA.

With the seed as written, guesses and slips come back within about 0.05 of the truth, except item 7’s slip (0.18 true, 0.27 estimated), on the item needing all three attributes. Each attribute is classified correctly for 94% to 95% of respondents, the whole profile for 84%. AIC prefers DINA (10,692) to G-DINA (10,712), as it should: DINA generated the data. Telling the model that item 4 needs only A1 (it needs A1 and A2) sends its slip from 0.19 to 0.50: respondents with A1 but not A2 now count as able, and get it wrong. The AIC rises to 10,828. Try n <- 300, or make some profiles rare in share.

Download this code to run it as large as you like in your own R session.

With real data

cdm_ecpe has two jobs: the latent class example, fit without a Q-matrix to ask whether its classes come out ordered, and the contrast for cognitive diagnosis, a test whose attributes may be one scale in disguise. frac20 is the main example, Tatsuoka’s fraction subtraction items. Both tables carry their Q-matrix in Qmatrix__ columns. Download the code; it needs poLCA, GDINA and mirt (Chalmers, 2012).

ECPE: classes without a Q-matrix

The Examination for the Certificate of Proficiency in English (ECPE) data are 28 grammar items answered by 2,922 respondents, from a tutorial by Templin and Hoffman (2013). The Q-matrix codes three attributes: morphosyntactic, cohesive and lexical rules (CDM package documentation; George et al., 2016). Keying: 1 is correct.

Code 
suppressPackageStartupMessages({library(poLCA); library(GDINA); library(mirt)})
options(digits = 7)  # R's default, in case a .Rprofile changes it
url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v60_0.cdm_ecpe/rows?format=csv"
ecpe <- read.csv(url)
# One row per response; the Q-matrix travels in the Qmatrix__ columns.
wide <- reshape(ecpe[, c("id", "item", "resp")], idvar = "id", timevar = "item",
                direction = "wide")
wide <- wide[order(wide$id), ]
e_items <- paste0("E", 1:28)
Xe <- as.matrix(wide[, paste0("resp.", e_items)])
colnames(Xe) <- e_items
Qe <- unique(ecpe[, c("item", "Qmatrix__skill1", "Qmatrix__skill2", "Qmatrix__skill3")])
Qe <- as.matrix(Qe[match(e_items, Qe$item), -1])
colnames(Qe) <- c("morphosyntactic", "cohesive", "lexical")
cat(nrow(Xe), "respondents,", ncol(Xe), "items,", sum(is.na(Xe)), "missing responses\n")
2922 respondents, 28 items, 0 missing responses
Code 
cat("Proportion correct: from", min(round(colMeans(Xe), 2)), "to", max(round(colMeans(Xe), 2)), "\n")
Proportion correct: from 0.43 to 0.9
Code 
# The same data ship with GDINA; check that the copies agree.
data(realdata_ECPE, package = "GDINA")
cat("Same responses as GDINA's copy:", all(Xe == as.matrix(realdata_ECPE$dat)),
    "  Same Q-matrix:", all(Qe == as.matrix(realdata_ECPE$Q)), "\n")
Same responses as GDINA's copy: TRUE   Same Q-matrix: TRUE

Nothing is missing, and proportions correct run from 0.43 to 0.90. First, set the attributes aside and fit one to five classes:

Code 
# poLCA wants categories coded 1, 2, ...; nrep = 5 starts each fit from five random
# places and keeps the best, since EM can stop on a local maximum.
d <- as.data.frame(Xe + 1)
f <- as.formula(paste("cbind(", paste(e_items, collapse = ","), ") ~ 1"))
set.seed(252)
lca <- lapply(1:5, function(k)
  poLCA(f, d, nclass = k, nrep = if (k == 1) 1 else 5, maxiter = 3000, verbose = FALSE))
rasch_e <- mirt(Xe, 1, "Rasch", verbose = FALSE)
twopl_e <- mirt(Xe, 1, "2PL", verbose = FALSE)
fits <- data.frame(
  model = c(paste(1:5, "classes"), "Rasch", "2PL"),
  logLik = round(c(sapply(lca, `[[`, "llik"), extract.mirt(rasch_e, "logLik"),
                   extract.mirt(twopl_e, "logLik")), 1),
  parameters = c(sapply(lca, `[[`, "npar"), extract.mirt(rasch_e, "nest"),
                 extract.mirt(twopl_e, "nest")),
  AIC = round(c(sapply(lca, `[[`, "aic"), extract.mirt(rasch_e, "AIC"),
                extract.mirt(twopl_e, "AIC"))),
  BIC = round(c(sapply(lca, `[[`, "bic"), extract.mirt(rasch_e, "BIC"),
                extract.mirt(twopl_e, "BIC"))))
print(fits, row.names = FALSE)
     model   logLik parameters   AIC   BIC
 1 classes -45376.9         28 90810 90977
 2 classes -42972.7         57 86059 86400
 3 classes -42547.8         86 85268 85782
 4 classes -42453.2        115 85136 85824
 5 classes -42402.6        144 85093 85954
     Rasch -42729.8         29 85518 85691
       2PL -42546.7         56 85205 85540

BIC chooses three classes (85,782); AIC keeps falling through five. But the 2PL has a lower BIC (85,540) than any class model, and with 56 parameters its log likelihood (−42,546.7) matches the three-class model’s (−42,547.8) with 86. Do the class profiles stack?

Code 
# Put each solution's classes in order of their average proportion correct, then
# count the items on which a higher class does worse than the class below it.
class_probs <- function(fit) {
  P <- sapply(fit$probs, function(p) p[, 2])   # classes x items: P(correct)
  o <- order(rowMeans(P))
  list(P = P[o, , drop = FALSE], share = fit$P[o])
}
for (k in 2:5) {
  cp <- class_probs(lca[[k]])
  out_of_order <- sum(apply(cp$P, 2, function(v) any(diff(v) < 0)))
  cat(sprintf("%d classes: shares %s; mean P(correct) %s; items out of order: %d\n", k,
              paste(sprintf("%.2f", cp$share), collapse = " "),
              paste(sprintf("%.2f", rowMeans(cp$P)), collapse = " "), out_of_order))
}
2 classes: shares 0.47 0.53; mean P(correct) 0.58 0.84; items out of order: 0
3 classes: shares 0.19 0.51 0.30; mean P(correct) 0.48 0.70 0.89; items out of order: 0
4 classes: shares 0.13 0.38 0.37 0.12; mean P(correct) 0.44 0.64 0.81 0.94; items out of order: 0
5 classes: shares 0.14 0.35 0.11 0.29 0.11; mean P(correct) 0.45 0.65 0.74 0.82 0.94; items out of order: 15

Up to four classes, yes, like the levels in the class-profiles widget. The five-class solution puts 15 items out of order. Where do its classes cross?

Code 
# The five-class solution, with its items grouped by the attributes the ECPE's
# Q-matrix says they need (1 = needed; order: morphosyntactic, cohesive, lexical).
cp5 <- class_probs(lca[[5]])
pattern <- apply(Qe, 1, paste, collapse = "")
tab5 <- aggregate(t(cp5$P), list(Q_row = pattern), mean)
names(tab5)[-1] <- paste0("class", 1:5)
tab5$items <- as.vector(table(pattern)[tab5$Q_row])
print(cbind(tab5[, 1, drop = FALSE], round(tab5[, 2:6], 2), items = tab5$items), row.names = FALSE)
 Q_row class1 class2 class3 class4 class5 items
   001   0.47   0.74   0.81   0.88   0.95    10
   010   0.57   0.74   0.66   0.87   0.93     4
   011   0.70   0.90   0.87   0.93   0.98     1
   100   0.42   0.52   0.68   0.71   0.93     5
   101   0.30   0.50   0.71   0.77   0.95     7
   110   0.60   0.76   0.77   0.89   0.97     1

Class 3 (11% of respondents) beats class 2 on items needing morphosyntactic rules (0.68 against 0.52 on the five needing only those, 0.71 against 0.50 on the seven also needing lexical rules) and trails it on the four needing only cohesive rules (0.66 against 0.74). That is a class that differs in which items it finds hard, and the difference follows the Q-matrix, which the fit never saw. But BIC doesn’t support a fifth class (85,954 against 85,782 for three). The unrestricted classes mostly come out ordered: a continuum cut into pieces.

Fraction subtraction: DINA

These are 536 students’ responses to the 20 items (Tatsuoka, 2002); the IRW table matches the copy in GDINA. Keying: 1 is correct.

Code 
url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v60_0.frac20/rows?format=csv"
frac <- read.csv(url)
wide <- reshape(frac[, c("id", "item", "resp")], idvar = "id", timevar = "item",
                direction = "wide")
wide <- wide[order(wide$id), ]
f_items <- paste0("item_", 1:20)
Xf <- as.matrix(wide[, paste0("resp.", f_items)])
colnames(Xf) <- f_items
Qf <- unique(frac[, c("item", paste0("Qmatrix__", 1:8))])
Qf <- as.matrix(Qf[match(f_items, Qf$item), -1])
colnames(Qf) <- paste0("A", 1:8)
cat(nrow(Xf), "respondents,", ncol(Xf), "items,", sum(is.na(Xf)), "missing responses\n")
536 respondents, 20 items, 0 missing responses
Code 
data(realdata_Tatsuoka1990, package = "GDINA")
cat("Same responses as GDINA's copy:", all(Xf == as.matrix(realdata_Tatsuoka1990$dat)),
    "  Same Q-matrix:", all(Qf == as.matrix(realdata_Tatsuoka1990$Q)), "\n")
Same responses as GDINA's copy: TRUE   Same Q-matrix: TRUE
Code 
cat("Attributes per item:", rowSums(Qf), "\n")
Attributes per item: 3 2 2 4 4 1 3 1 1 4 3 2 4 2 2 2 3 4 5 4
Code 
cat("Items needing each attribute:", colSums(Qf), "\n")
Items needing each attribute: 3 13 3 5 8 2 19 3

Items need one to five attributes. A7 is needed by 19 items and A2 by 13; A1, A3 and A8 by three each, and A6 by two. Now DINA:

Code 
dina_f <- GDINA(Xf, Qf, model = "DINA", verbose = 0)
gs <- coef(dina_f, "gs")
print(round(t(gs[, c("guessing", "slip")]), 2))
         Item 1 Item 2 Item 3 Item 4 Item 5 Item 6 Item 7 Item 8 Item 9 Item 10
guessing   0.03   0.02   0.00   0.22   0.30   0.10   0.03   0.44   0.29    0.03
slip       0.09   0.04   0.13   0.11   0.17   0.04   0.20   0.18   0.25    0.21
         Item 11 Item 12 Item 13 Item 14 Item 15 Item 16 Item 17 Item 18
guessing    0.07    0.13    0.01    0.06    0.03    0.11    0.04    0.12
slip        0.08    0.04    0.33    0.06    0.11    0.11    0.14    0.14
         Item 19 Item 20
guessing    0.02    0.01
slip        0.24    0.16
Code 
cat(sprintf("Items with slip at most 0.25: %d of 20. Largest guess: %.2f (item %d)\n",
            sum(gs[, "slip"] <= 0.25), max(gs[, "guessing"]), which.max(gs[, "guessing"])))
Items with slip at most 0.25: 19 of 20. Largest guess: 0.44 (item 8)
Code 
cat("Log likelihood", round(logLik(dina_f), 1), "with", npar(dina_f)$`No. of parameters`,
    "parameters\n")
Log likelihood -4402.3 with 295 parameters

Slips are at most 0.25 on 19 items; item 13’s is 0.33. Most guesses are near zero, as expected for free-response subtraction. Item 8 stands out with a guess of 0.44: it needs only A7, yet students without A7 answer it correctly 44% of the time by the model’s estimate. Perhaps it can be solved another way, or its Q row misses something. Who has which attributes?

Code 
prev <- extract(dina_f, "prevalence")$all[, "Level1"]
cat("Proportion mastering each attribute:", round(prev, 2), "\n")
Proportion mastering each attribute: 0.58 0.77 0.72 0.69 0.6 0.79 0.81 0.82
Code 
att_f <- personparm(dina_f, "EAP")   # attribute mastered if its posterior probability > 0.5
n_att <- rowSums(att_f)
cat("Respondents with all eight:", sum(n_att == 8), sprintf("(%.1f%%)", 100 * mean(n_att == 8)),
    "  Distinct profiles among", nrow(att_f), "respondents:",
    length(unique(apply(att_f, 1, paste, collapse = ""))), "\n")
Respondents with all eight: 200 (37.3%)   Distinct profiles among 536 respondents: 49

Between 58% (A1) and 82% (A8) of students master each attribute, and 200 of 536 (37.3%) are classified as having all eight. Counting an attribute as mastered when its posterior probability is above 0.5, the students fall into 49 profiles.

Recall from The Rasch model: under the Rasch model the sum score is sufficient for \(\theta\). Two respondents with the same number right get the same estimate, whichever items they got right.

DINA reads which items. Among students with the same sum score, do the profiles differ?

Code 
r <- rowSums(Xf)
profile <- apply(att_f, 1, paste, collapse = "")
by_r <- data.frame(sum_score = sort(unique(r)),
                   respondents = as.vector(table(r)),
                   profiles = as.vector(tapply(profile, r, function(p) length(unique(p)))),
                   fewest_attributes = as.vector(tapply(n_att, r, min)),
                   most_attributes = as.vector(tapply(n_att, r, max)))
print(by_r[by_r$sum_score %in% c(0, 5, 10, 15, 20), ], row.names = FALSE)
 sum_score respondents profiles fewest_attributes most_attributes
         0          13        1                 0               0
         5          31       11                 2               5
        10          16       10                 5               7
        15          23        3                 7               8
        20          30        1                 8               8
Code 
cat("\nThe 16 respondents with 10 correct, by profile (A1 ... A8):\n")

The 16 respondents with 10 correct, by profile (A1 ... A8):
Code 
print(sort(table(profile[r == 10]), decreasing = TRUE))

01010111 11010111 01110111 00111111 01011110 01101011 01110011 10010111 
       4        3        2        1        1        1        1        1 
11101111 11110111 
       1        1
Code 
rasch_f <- mirt(Xf, 1, "Rasch", verbose = FALSE)
theta_f <- fscores(rasch_f)[, 1]
cat("\nCorrelation of attributes mastered with the sum score:", round(cor(n_att, r), 2),
    "; with Rasch theta:", round(cor(n_att, theta_f), 2), "\n")

Correlation of attributes mastered with the sum score: 0.95 ; with Rasch theta: 0.93

At a sum score of 10, 16 students fall into 10 different profiles, with five to seven attributes; at 0 and at 20 everyone shares a profile. The number of attributes mastered correlates 0.95 with the sum score and 0.93 with the Rasch \(\theta\). The count of attributes mostly says what the sum score says; the profile adds which skills lie behind it, and that addition rests on the Q-matrix. How do the models compare?

Code 
gdina_f <- GDINA(Xf, Qf, model = "GDINA", verbose = 0)
ho_f <- GDINA(Xf, Qf, model = "DINA", att.dist = "higher.order", verbose = 0)
twopl_f <- mirt(Xf, 1, "2PL", verbose = FALSE)
cdm_row <- function(name, m) data.frame(model = name, logLik = round(as.numeric(logLik(m)), 1),
  parameters = npar(m)$`No. of parameters`, AIC = round(AIC(m)), BIC = round(BIC(m)))
irt_row <- function(name, m) data.frame(model = name, logLik = round(extract.mirt(m, "logLik"), 1),
  parameters = extract.mirt(m, "nest"), AIC = round(extract.mirt(m, "AIC")),
  BIC = round(extract.mirt(m, "BIC")))
print(rbind(cdm_row("DINA", dina_f), cdm_row("G-DINA", gdina_f),
            cdm_row("higher-order DINA", ho_f), irt_row("Rasch", rasch_f),
            irt_row("2PL", twopl_f)), row.names = FALSE)
             model  logLik parameters  AIC   BIC
              DINA -4402.3        295 9395 10658
            G-DINA -4265.8        445 9422 11328
 higher-order DINA -4509.6         48 9115  9321
             Rasch -4796.3         21 9635  9725
               2PL -4640.1         40 9360  9532
Code 
# The higher-order model: each attribute is a Rasch-type item on one theta (slope
# fixed at 1, GDINA's default), so its profile distribution has 8 parameters, not 255.
cat("\nHigher-order DINA, attribute intercepts:", round(coef(ho_f, "lambda")[, "intercept"], 2),
    "; slopes:", unique(coef(ho_f, "lambda")[, "slope"]), "\n")

Higher-order DINA, attribute intercepts: 0.44 2.48 2.86 1.21 0.58 3.6 1.39 2.54 ; slopes: 1
Code 
# A Q-matrix with no thought behind it: redraw the entries for attributes 5 to 8 at
# random, keeping each attribute's share of items (as in 252's c7/3_cdm.R).
set.seed(10103101)
Q_random <- Qf
for (k in 5:8) Q_random[, k] <- rbinom(20, 1, mean(Qf[, k]))
dina_random <- GDINA(Xf, Q_random, model = "DINA", verbose = 0)
cat("\nDINA with attributes 5-8 of the Q-matrix redrawn at random: AIC", round(AIC(dina_random)), "\n")

DINA with attributes 5-8 of the Q-matrix redrawn at random: AIC 9769

DINA beats the Rasch model by AIC (9,395 against 9,635) but not the 2PL (9,360), and by BIC trails both: 255 of its 295 parameters are class shares. G-DINA’s extra parameters don’t pay (AIC 9,422). The higher-order DINA wins on both criteria (AIC 9,115, BIC 9,321, 48 parameters), with all its attributes riding on one \(\theta\). The Q-matrix carries information: with attributes 5 to 8 redrawn at random, the AIC rises to 9,769. Do the data suggest a better one?

Code 
# Empirical Q-matrix validation by the proportion of variance accounted for (PVAF;
# de la Torre & Chiu, 2016), starting from the DINA fit and from the G-DINA fit.
# For each suggestion, refit the same model with the suggested Q-matrix.
for (m in list(list("DINA", dina_f), list("GDINA", gdina_f))) {
  Q_sug <- extract(Qval(m[[2]], method = "PVAF", eps = 0.95), "sug.Q")
  changed <- which(rowSums(Q_sug != Qf) > 0)
  refit <- GDINA(Xf, Q_sug, model = m[[1]], verbose = 0)
  cat(sprintf("From %s: changes suggested to %d items (%s). AIC %.0f with the original Q, %.0f with the suggested one.\n",
              m[[1]], length(changed), paste(changed, collapse = ", "), AIC(m[[2]]), AIC(refit)))
}
From DINA: changes suggested to 10 items (4, 5, 7, 8, 9, 10, 13, 16, 19, 20). AIC 9395 with the original Q, 9527 with the suggested one.
From GDINA: changes suggested to 8 items (1, 4, 7, 8, 9, 12, 13, 19). AIC 9422 with the original Q, 9340 with the suggested one.

Starting from the DINA fit, PVAF suggests changing 10 of the 20 items, and DINA refit with the suggested matrix is worse (AIC 9,527 against 9,395). Starting from G-DINA, the model the method is defined under, it suggests 8, and G-DINA with the suggested matrix improves (9,340 against 9,422). The suggestions depend on the starting model; a refit and a look at the item are what test each one (problem 3).

ECPE with its Q-matrix

Back to the ECPE, now with its three attributes and eight profiles.

Code 
gdina_e <- GDINA(Xe, Qe, model = "GDINA", verbose = 0)
shares <- extract(gdina_e, "posterior.prob")[1, ]
print(round(sort(shares, decreasing = TRUE), 3))
  111   000   011   001   101   110   010   100 
0.349 0.302 0.180 0.125 0.016 0.014 0.012 0.002
Code 
cat("Four largest profiles together:", round(sum(sort(shares, decreasing = TRUE)[1:4]), 3), "\n")
Four largest profiles together: 0.956

Four profiles hold 95.6% of respondents, the other four under 2% each. The four form a chain: lexical rules first, then cohesive, then morphosyntactic, so knowing how many attributes a respondent has tells you which. Templin and Bradshaw (2014) used these data to test for such a hierarchy. With only the four profiles allowed:

Code 
dina_e <- GDINA(Xe, Qe, model = "DINA", verbose = 0)
# G-DINA allowed only the four profiles of a chain: none, lexical, lexical + cohesive, all
chain <- matrix(c(0, 0, 0,  0, 0, 1,  0, 1, 1,  1, 1, 1), ncol = 3, byrow = TRUE)
gdina_chain <- GDINA(Xe, Qe, model = "GDINA", att.str = chain, verbose = 0)
print(rbind(cdm_row("DINA", dina_e), cdm_row("G-DINA", gdina_e),
            cdm_row("G-DINA, four-profile chain", gdina_chain),
            irt_row("Rasch", rasch_e), irt_row("2PL", twopl_e)), row.names = FALSE)
                      model   logLik parameters   AIC   BIC
                       DINA -42841.6         63 85809 86186
                     G-DINA -42738.6         81 85639 86124
 G-DINA, four-profile chain -42751.3         68 85639 86045
                      Rasch -42729.8         29 85518 85691
                        2PL -42546.7         56 85205 85540

The chain matches the full G-DINA’s AIC (85,639) with 13 fewer parameters, so its BIC is better (86,045 against 86,124). But the Rasch model (AIC 85,518) and the 2PL (85,205) beat every CDM on both criteria. A four-step chain is an ordinal scale, and here the continuum describes the responses better. Both routes through the ECPE, classes without a Q-matrix and profiles with one, arrive at ordered groups. In the fraction data the Q-matrix adds information, but the model that uses it best ties the attributes to one continuum.

Whether a test’s attributes behave as separate skills or collapse toward one continuum, as in both datasets here, is hard to know in advance, and AIC and BIC can disagree about it. Responses the models haven’t seen are the acid test: fit a CDM and a unidimensional model, compare how well each predicts held-out responses (problem 3), and let that comparison guide the choice.

Problems

  1. The DINA likelihood. Write the likelihood of one respondent’s responses under DINA as a latent class model over the \(2^K\) profiles. Show that with a single attribute (\(K = 1\)), DINA is a two-class latent class model whose item probabilities are \(g_i\) in one class and \(1 - s_i\) in the other. How many parameters does each model have, and are they the same model?
  2. The Q-matrix in an LLTM. Fit resp ~ 0 + Qmatrix__1 + ... + Qmatrix__8 + (1 | id) to frac20 with glmer, as in What is an item? (PS7#2 B). Do items with the same Q row have similar Rasch difficulties? What does the LLTM assume about attributes that DINA doesn’t, and the reverse?
  3. Held out. Hold out a random 20% of frac20 responses (set them to missing), fit the Rasch model, the 2PL, DINA and the higher-order DINA to the rest, and compute each model’s IMV over the Rasch model on the held-out responses (Model fit and out-of-sample prediction). Does the ranking agree with AIC, and with what Ma et al. (2026) found for these data? The IRW vignette Calculating IMV values has the cross-validation code. Then run the full PVAF analysis: which suggested changes make sense for the items as written? (PS7#2 D.)
  4. Skill reports from a state test. A district wants attribute-level reports for each student from a state mathematics test built from a blueprint of content strands, not a Q-matrix. Advise them. What would you want to see before any report goes home?
  5. Design a Q-matrix. Write a Q-matrix for eight items on three attributes in which each attribute has an item that needs it alone. Why does that matter for telling the \(2^3\) profiles apart under DINA? Try a Q-matrix without such items in the simulation.
  6. Challenge (open). When attributes form a hierarchy, as the ECPE’s appear to, is a CDM a new kind of model or a unidimensional model with cut points? Set out what each reading predicts about data you could collect, such as items needing a later attribute but not an earlier one. What would convince you?

Stuck, or curious about something this lesson doesn’t cover? Ask Claude about this lesson opens a new conversation that already knows what the lesson is about; add your question at the end. Claude shows a caution notice on any prompt that arrives through a link; that’s expected.

Going further

  • Rupp, A. A., Templin, J., & Henson, R. A. (2010). Diagnostic measurement: Theory, methods, and applications. Guilford Press. https://openlibrary.org/works/OL15026394W. A book-length treatment of cognitive diagnosis models, from Q-matrices to estimation.
  • Ma, W., & de la Torre, J. (2020). GDINA: An R package for cognitive diagnosis modeling. Journal of Statistical Software, 93(14). https://doi.org/10.18637/jss.v093.i14. The package used here, with examples.
  • Linzer, D. A., & Lewis, J. B. (2011). poLCA: An R package for polytomous variable latent class analysis. Journal of Statistical Software, 42(10). https://doi.org/10.18637/jss.v042.i10. Latent class models in R, including covariates for class membership.
  • Lindsay, B., Clogg, C. C., & Grego, J. (1991). Semiparametric estimation in the Rasch model and related exponential response models, including a simple latent class model for item analysis. Journal of the American Statistical Association, 86(413), 96–107. https://doi.org/10.1080/01621459.1991.10475008. Why a few ordered latent classes can stand in for the Rasch model’s continuum.
  • de la Torre, J., & Douglas, J. A. (2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69(3), 333–353. https://doi.org/10.1007/bf02295640. Attributes driven by one continuous trait.
  • DeCarlo, L. T. (2011). On the analysis of fraction subtraction data: The DINA model, classification, latent class sizes, and the Q-matrix. Applied Psychological Measurement, 35(1), 8–26. https://doi.org/10.1177/0146621610377081. DINA fit as a latent class model to the fraction data, and how it classifies students who get every item wrong.
  • Templin, J., & Bradshaw, L. (2014). Hierarchical diagnostic classification models: A family of models for estimating and testing attribute hierarchies. Psychometrika, 79(2), 317–339. https://doi.org/10.1007/s11336-013-9362-0. Testing an attribute hierarchy, with a grammar test as the example.
  • Templin, J., & Bradshaw, L. (2013). Measuring the reliability of diagnostic classification model examinee estimates. Journal of Classification, 30(2), 251–275. https://doi.org/10.1007/s00357-013-9129-4. How consistently a CDM classifies respondents, a question this lesson leaves open.
  • Zhou, W., Liu, Y., Tang, R., & Wilson, M. (2026). Aligning LLM-simulated and human examinees for psychometric calibration: A cognitive diagnostic profiling approach. arXiv. https://arxiv.org/abs/2607.26317. Attribute profiles written as student descriptions, used to prompt language models to answer the fraction items (a preprint).
  • Ma, W. A., Liu, Y., Kanopka, K., Ma, W., & Domingue, B. W. (2026). A comparison of the predictive performance of continuous and class-based latent trait models. Psychological Methods. Advance online publication. https://doi.org/10.1037/met0000860. CDMs against the 2PL on held-out responses, in simulations and 11 datasets.
  • The InterModel Vigorish: IRT versus cognitive diagnostic models, an interactive simulation of when a CDM predicts better than the 2PL.
  • Rost, J. (1990). Rasch models in latent classes: An integration of two approaches to item analysis. Applied Psychological Measurement, 14(3), 271–282. https://doi.org/10.1177/014662169001400305. Mixture IRT: a continuum within each class.
  • von Davier, M., & Rost, J. (2016). Logistic mixture-distribution response models. In W. J. van der Linden (Ed.), Handbook of item response theory: Vol. 1. Models (pp. 393–406). Chapman and Hall/CRC. https://doi.org/10.1201/9781315374512. A handbook chapter on mixture item response models.
  • Wang, C., & Chang, H.-H. (2017). Cognitive diagnostic assessment. In W. J. van der Linden (Ed.), Handbook of item response theory: Vol. 3. Applications (pp. 329–348). Chapman and Hall/CRC. https://doi.org/10.1201/9781315117430. A handbook chapter on cognitive diagnostic assessment.

Data sources

Data from the Item Response Warehouse (each table name links to its IRW page); references from IRW biblio.

  • frac20: Tatsuoka, C. (2002). Data analytic methods for latent partially ordered classification models. Journal of the Royal Statistical Society: Series C (Applied Statistics), 51(3), 337-350. https://doi.org/10.1111/1467-9876.00272. Licence: GPL-3.0.
  • cdm_ecpe: Templin, J., & Hoffman, L. (2013). Obtaining diagnostic classification model estimates using Mplus. Educational Measurement: Issues and Practice, 32, 37-50. https://doi.org/10.1111/emip.12010. Licence: GPL-3.0.

For instructors

  • Source: EDUC 252 class 7 (the c7 slides are missing from the export, so the lesson is built from the code) and problem set 7 (#2), with the code in ben-domingue/252: c7/3_cdm.R and ps7/cdm.R.
  • Timing: the core ideas and widgets take about 50 minutes; the real-data section is a 30-minute lab. The latent class fits take a minute or two, so run them before class.
  • Changed from 252: latent class analysis now comes first, as the general model; the ECPE’s unrestricted classes and its four-profile chain are new. The DINA fit to the fraction data follows c7/3_cdm.R (same data and Q-matrix, and the same log likelihood, −4,402.3), including its randomly redrawn Q-matrix. Attributes are classified by posterior probability above 0.5 rather than by the most probable profile, because personparm(..., "MAP") breaks ties at random.
  • Held back: longitudinal CDMs (the IRW’s hmcdm_spatialreasoning table), mixture IRT (in Going further) and the classification accuracy of attribute profiles are not taken up in this course. The held-out comparison and the full PVAF analysis are problem 3. cdm_timss07 (15 attributes) is too slow for class.
  • Verdict: whether the attributes are separate skills or one continuum is hard to know in advance; compare a CDM and a unidimensional model on held-out responses and let that guide the choice (F27, Ben 09-30).
  • Solutions: held back in solutions/cdm.qmd (not published). Problem 6 is open.