The Dutch Identity: Manifest Probabilities and IRT Models
A tour of Holland’s (1990) Dutch Identity — the result that log-manifest probabilities of dichotomous response patterns are approximately log-quadratic in the item responses. We ground each theoretical claim in IRW data: enumerating manifest probabilities directly, fitting loglinear models, diagnosing the Rasch model via total-score frequencies, assessing dimensionality from second-order interactions, and testing Holland’s conservation law for model complexity.
Published
October 23, 2026
Note
This vignette is currently in development. Content and results may change.
Background
Holland (1990) introduced the Dutch Identity — a result about the mathematical structure that any unidimensional IRT model imposes on the joint distribution of dichotomous item responses. The key insight is that the manifest (unconditional) probability of a response pattern \(\mathbf{x} =
(x_1, \ldots, x_J) \in \{0,1\}^J\) can be written as an exponential family whose sufficient statistics are functions of the item responses. For large tests with smooth item response functions (IRFs), this reduces to a strikingly simple approximation:
This is log-quadratic in the \(x_j\), requires only \(2J + 1\) parameters regardless of \(J\), and makes precise predictions about the structure of second-order interactions between items. The sections below work through the implications one step at a time, using IRW datasets to make each idea concrete.
Note
All figures use pre-computed results from dutch_identity_compute.R. The dataset used in Sections 1–3 is ALSECYPIAMH_WU_2022_NEI (\(J = 5\) items, \(N = 7,842\) respondents, computed October 2026).
Section 1: The manifest probabilities directly
For a test of \(J\) binary items, a respondent’s answers define a response pattern \(\mathbf{x} \in \{0,1\}^J\). There are \(2^J\) possible patterns, each occurring with manifest probability \(p(\mathbf{x})\). Before modeling, it is worth asking: what do these probabilities actually look like? Most patterns are extremely rare; a handful dominate. Holland (1990) motivates the loglinear approach precisely because the \(p(\mathbf{x})\) are too numerous and sparse to model directly for large \(J\).
For \(J \leq 10\) items and \(N \geq 2{,}000\) respondents the full distribution is enumerable. Below we tabulate every \(2^J\) pattern and display the empirical \(p(\mathbf{x})\).
Figure 1: Distribution of empirical manifest probabilities across all 2^J response patterns (log scale). Most patterns are never observed (grey); a few common patterns carry most of the probability mass. Hover over a point to see its response pattern.
Code
patterns |>filter(observed) |>ggplot(aes(x =factor(total), y = p_obs, fill = p_obs)) +geom_col(width =0.7, show.legend =FALSE) +scale_fill_gradient(low ="#d1e5f0", high = irw_blue) +labs(x ="Total score k",y ="p(x)",caption ="Each bar is the total probability mass at that score level." )
Figure 2: Total probability mass at each sum score k = Σ xⱼ. Patterns within each score group contribute additively.
What to notice. Probability mass is extremely concentrated: the most common pattern can account for a large share of the total. The long tail of near-zero patterns means a \(2^J\)-parameter saturated model is hopeless — most cells have zero or one observation. This motivates the loglinear parameterisation in Section 2.
The bits-per-item vignette uses this same distribution to count how much a test tells you about a person: the entropy of p(x), minus the part that is chance rather than differences between people. It also shows that the independence model of Section 2 falls short by exactly the information that items share.
Section 2: Log-manifest probabilities and the independence baseline
Under item independence — no latent variable, no associations — \(\log
p(\mathbf{x})\) is a linear function of the item responses:
This independence loglinear model has \(J + 1\) parameters. It is equivalent to a Poisson regression of pattern counts on the \(J\) binary item predictors. When it fits, all structure in the data is captured by item marginals alone.
Any IRT model implies that \(\log p(\mathbf{x})\)deviates from this baseline — the latent variable induces positive associations among items that the main-effects model misses. Plotting observed versus fitted \(\log p(\mathbf{x})\) makes those deviations visible.
Figure 3: Observed vs. fitted log p(x) under the independence loglinear model. Points on the diagonal are fully explained by item marginals; points below the line are more probable than independence predicts (positive local dependence).
What to notice. Systematic curvature away from the diagonal — especially for extreme total-score patterns — signals that main effects alone are not enough. The points that deviate most are those where all or nearly all items are answered the same way: item correlation inflates these relative to the independence prediction. This is exactly the second-order structure that Holland (1990) formalises.
Section 3: Fitting the second-order exponential model
The full second-order exponential (SOE) model augments the independence model with all \(\binom{J}{2}\) pairwise interaction terms:
This is still a Poisson GLM, now with \(1 + J + \binom{J}{2}\) parameters. The \(J \times J\) matrix \(\boldsymbol{\Gamma}\) (with \(\Gamma_{jk} = \gamma_{jk}\) and \(\Gamma_{jj} = 0\)) encodes all pairwise item associations after removing main effects.
The Dutch Identity (Holland 1990, Theorem 1) predicts that \(\boldsymbol{\Gamma}\) is approximately rank-1: \(\Gamma_{jk} \approx \gamma_j
\gamma_k\). Under this constraint, the interaction structure collapses to:
Table 1: Log-likelihood and BIC for three nested loglinear models, all fitted by Poisson maximum likelihood to the same table of 2^J pattern counts with zero cells included. The restricted log-quadratic uses 2J + 1 parameters; the full SOE uses 1 + J + J(J-1)/2. BIC uses n = N respondents for every row.
Model
Parameters
Log-lik
BIC
ΔBIC vs. best
Independence
6
-3816.7
7687.1
7260.5
Log-quadratic (rank-1 Γ)
11
-199.3
497.3
70.7
Full SOE
16
-141.6
426.6
0.0
All three models are fitted to the same \(2^J\) pattern counts — zero cells included — by Poisson maximum likelihood, so the log-likelihoods are on one scale, and BIC uses \(n = N = 7,842\) respondents for every row. Each model has a free intercept, so its fitted counts sum to \(N\); the multinomial log-likelihoods therefore differ from the Poisson ones by a constant that is the same for all three (-24523.9, -20906.6, -20848.8 for independence, log-quadratic and SOE respectively).
The three models are nested: independence is the log-quadratic with \(\gamma = 0\), and the log-quadratic is the full SOE with \(\Gamma_{jk}\) forced to \(\gamma_j \gamma_k\). The log-likelihoods must therefore increase down the table, and they do. The rank-1 fit is by direct Poisson maximum likelihood rather than by a GLM, since \(\gamma\) enters only through \((\mathbf{x}'\gamma)^2\) and so is identified only up to sign; of 20 starting values, 16 reached the reported optimum.
The key question is whether the restricted log-quadratic — which enforces rank-1 structure on \(\boldsymbol{\Gamma}\) — captures most of the association signal while using only \(1 + 2J\) parameters instead of \(1 + J + \binom{J}{2}\). It captures nearly all of it: the rank-1 model recovers 98.4% of the independence-to-SOE log-likelihood gap, \((\ell_{\text{rq}} - \ell_{\text{ind}})/(\ell_{\text{soe}} - \ell_{\text{ind}})\), with 5 fewer parameters. BIC still prefers the full SOE, but only by 70.7 — small beside the 7260.5 that separates independence from the SOE.
Figure 4: Observed vs. fitted log-probability for each of the three loglinear models across all observed response patterns. Points on the dashed diagonal indicate perfect fit; scatter above the diagonal means the model over-predicts that pattern.
Structure of the Γ matrix
Code
gmat <- s1_3$gamma_matrixgmat_long <-as.data.frame(as.table(gmat)) |>rename(item1 = Var1, item2 = Var2, gamma = Freq) |>filter(item1 != item2)ggplot(gmat_long, aes(x = item2, y = item1, fill = gamma)) +geom_tile(colour ="white") +scale_fill_gradient2(low = irw_red, mid ="white", high = irw_blue,midpoint =0, name ="γⱼₖ" ) +labs(x =NULL, y =NULL) +theme(axis.text.x =element_text(angle =45, hjust =1),panel.grid =element_blank() )
Figure 5: Heatmap of the J×J interaction matrix Γ from the full SOE model. Under the Dutch Identity, off-diagonal entries should follow a rank-1 pattern γⱼγₖ — all rows should look like rescaled versions of each other.
Under a rank-1 \(\boldsymbol{\Gamma}\), every row is a scalar multiple of every other — the color profile across columns should look identical up to intensity. Systematic blocks or rows that deviate from the common pattern signal local dependence that a single latent factor cannot absorb.
Is Γ rank-1?
The Dutch Identity’s prediction about \(\boldsymbol{\Gamma}\) concerns the off-diagonal entries only. Expanding the quadratic term and using \(x_j^2 = x_j\) for binary items,
so the squared terms are absorbed into the main effects and the prediction is \(\Gamma_{jk} \approx \gamma_j \gamma_k\) for \(j \neq k\). The diagonal is free; in the SOE fit \(\Gamma_{jj} = 0\) by construction, because there is no \(x_j^2\) term to estimate.
That construction has a consequence worth stating before any eigenvalues are plotted: \(\operatorname{trace}(\boldsymbol{\Gamma}) = 0\), so the eigenvalues of \(\boldsymbol{\Gamma}\) must sum to zero. Large negative eigenvalues are forced by the zero diagonal and say nothing about dimensionality. Indeed for any\(\boldsymbol{\gamma}\), the exactly rank-1 structure \(\boldsymbol{\gamma}\boldsymbol{\gamma}' - \operatorname{diag}(\gamma_j^2)\) has at most one positive eigenvalue, from which \(\lambda_1 / \sum_j |\lambda_j| \le 1/2\) with equality precisely when one eigenvalue is positive. A scree plot of \(\boldsymbol{\Gamma}\) is therefore not a test of the rank-1 prediction, and “the leading eigenvalue explains 50% of the absolute eigenvalue mass” is the ceiling of that statistic rather than a shortfall. An earlier version of this page read it the other way around.
The diagnostic that does test the prediction leaves the diagonal alone: fit the best rank-1 structure to the off-diagonal entries, minimising \(\sum_{j<k} w_{jk}(\Gamma_{jk} - g_j g_k)^2\) with weights \(w_{jk} = 1/\mathrm{SE}(\Gamma_{jk})^2\) from the SOE fit, and ask what share of the off-diagonal (weighted) sum of squares it explains.
Figure 6: Off-diagonal entries of Γ from the full SOE fit against γⱼγₖ from the fitted restricted log-quadratic model, with the 45° line. Error bars are ±1 SE of the SOE coefficient. This is the direct check of the Dutch Identity’s prediction; the diagonal plays no part in it.
Code
eig_df <-tibble(component =seq_along(s1_3$eigenvalues),eigenvalue = s1_3$eigenvalues,reference = s1_3$gamma_rank1$reference_eigen)ggplot(eig_df, aes(x = component, y = eigenvalue)) +geom_col(fill = irw_blue, alpha =0.7, width =0.6) +geom_hline(yintercept =0, colour = irw_grey) +geom_point(aes(y = reference), colour = irw_red, size =2.6) +scale_x_continuous(breaks =seq_len(s1_3$J)) +labs(x ="Component",y ="Eigenvalue of Γ",caption =paste0("Bars: observed Γ. Points: exactly rank-1 reference. Both sum to zero ","by construction (trace Γ = 0)." ) )
Figure 7: Eigenvalues of Γ (bars), with the eigenvalues of an exactly rank-1 reference structure γγ′ with its diagonal zeroed (points). Both sets are constrained to sum to zero by the zero diagonal, so the negative eigenvalues are not evidence against rank-1 — the rank-1 ideal has them too.
What to notice. The restricted log-quadratic model (\(1 + 2J\) parameters) closes 98.4% of the gap between the independence model and the full SOE (\(1 + J + \binom{J}{2}\) parameters) — the Dutch Identity’s economy claim, holding up here.
The rank-1 claim about \(\boldsymbol{\Gamma}\) itself is a separate question, and it needs a reference point: a test this short at this \(N\) will never return an exactly rank-1 \(\boldsymbol{\Gamma}\), so a high share proves nothing on its own. We simulate 200 sets of pattern counts from the fitted restricted log-quadratic model at the same \(N\) — a world where the Dutch Identity holds exactly — refit the SOE, and recompute the diagnostic. Under that null the weighted off-diagonal share has median 99.7% and 95% of replicates fall in [99.3%, 100.0%].
The observed share is 94.1% (94.0% unweighted), below every one of the 200 null replicates. So \(\boldsymbol{\Gamma}\) here is not rank-1, and the departure is not a sampling artefact. Figure 6 shows where it lives: 5 of the 10 off-diagonal entries sit more than three standard errors from \(\gamma_j\gamma_k\), the largest being item_NEI2–item_NEI5 at -6.4 SE. That is local dependence between particular item pairs — structure a single latent dimension cannot absorb.
Note also that the eigenvalue count has no power to detect this. \(\boldsymbol{\Gamma}\) has exactly 1 positive eigenvalue, which is the rank-1 prediction — but so do 100% of the null replicates, and so does any rank-1 matrix with a zeroed diagonal. Passing that check is not evidence.
The two findings are not in conflict, and the combination is the useful lesson. Rank-1 is rejected as a description of \(\boldsymbol{\Gamma}\) while remaining an excellent predictive approximation: the five extra interaction parameters the SOE spends on the non-rank-1 structure buy it only 1.6% of the likelihood gap. Statistical rejection and practical adequacy are answers to different questions, and at \(N = 7,842\) a real but small departure is detectable without being important for prediction.
Section 4: The γₖ diagnostic for the Rasch model
Holland (1990) Corollary 2 gives a striking prediction about the Rasch model. If ability \(\theta \sim N(\mu, \sigma^2)\), the probability of obtaining total score \(k\) is:
where \(\bar{y}_+ = \sum_j \bar{p}_j\) is the expected total score at the mean ability. In other words, \(\log \pi_k\) is a quadratic function of \(k\), centred on \(\bar{y}_+\). Systematic deviations indicate either non-normality of the ability distribution or Rasch misfit.
We compare two IRW datasets, chosen on the quantity these figures actually show — the \(R^2\) of the quadratic fit to \(\log \pi_k\). One where that quadratic holds closely (cdm_ecpe, \(J = 28\), \(N = 2,922\), \(R^2 = 0.988\)) and one where it does not (gilbert_meta_100, \(J = 12\), \(N = 2,118\), \(R^2 = 0.56\)).
Rasch misfit itself is quoted per response, \(\Delta\ell / (NJ)\) for the 2PL against the Rasch model: 0.0022 and 0.0043 respectively. Raw \(\Delta\ell\) is not comparable across datasets — it grows with the number of responses, so a larger dataset looks like worse Rasch fit for free, and the raw values here (183.1 and 110) would overstate the contrast by a wide margin.
Code
plot_gk(s4$good_rasch$gamma_k_data)
Figure 8: Total-score log-frequencies for the dataset with the closest quadratic fit. Points follow a smooth quadratic, consistent with Holland’s Corollary 2.
The quadratic \(R^2\) is 0.988 for this dataset.
Code
plot_gk(s4$poor_rasch$gamma_k_data)
Figure 9: Total-score log-frequencies for the dataset with the worst quadratic fit. Deviations from the quadratic suggest a non-normal ability distribution or non-constant discrimination.
The quadratic \(R^2\) is 0.56 for this dataset.
What to notice. When Rasch holds, \(\log \pi_k\) traces a clean inverted parabola. Deviations can take two forms: (1) a sharp peak at the mode with thin tails, suggesting a more leptokurtic ability distribution; or (2) asymmetry or a kink, consistent with varying item discriminations. The \(R^2\) from the quadratic regression is a quick, model-free diagnostic, but it is a comparative one rather than a test against a fixed cutoff. Across the 8 datasets screened here it ranges from 0.56 to 0.99; the top of that range is about as quadratic as real total-score distributions get, and values below roughly 0.95 are worth investigating.
The two diagnostics need not agree, and here they do not. Corollary 2 is a joint prediction about Rasch fit and normality of the ability distribution, so either can break it. cdm_ecpe happens to be best on both counts — highest quadratic \(R^2\) and the smallest per-response 2PL gain in the set. But gilbert_meta_100, with much the worst quadratic fit, is only middling on Rasch misfit per response (0.0043, against a worst-in-set 0.021 for anunciacao_2024_intelligence_gmi). What its \(\log \pi_k\) curve is flagging is the shape of the score distribution, not non-constant discrimination — which is the caveat in the box at the foot of this page, met in the data.
Section 5: Dimensionality via second-order interaction matrices
Holland ((1990), Section 3.3) shows that for a \(D\)-dimensional IRT model, the off-diagonal structure of \(\boldsymbol{\Gamma}\) is approximately of rank \(D\) — and, as in Section 3, this is an approximation that sharpens as \(J\) grows, about the off-diagonals with the diagonal left free. It is tempting to read it as “count the non-negligible eigenvalues of \(\boldsymbol{\Gamma}\)”, but that reading does not work: the SOE fit forces \(\Gamma_{jj} = 0\), so the eigenvalues sum to zero and the negative ones are an artefact of the constraint rather than a count of dimensions. Section 3 gives the diagnostic that does work, and shows that even the \(D = 1\) case has to be judged against a simulated rank-1 reference rather than by eye. So this section stays with the standard approach — eigenvalues of the inter-item correlation matrix, whose diagonal is 1 and which is a different, conventional setup — and treats \(\boldsymbol{\Gamma}\) only through Section 3’s machinery.
We use ALSECYPIAMH_WU_2022_NEI (\(J = 5\)) as the unidimensional candidate and CPDMMC_Kunnari_2020_PDP (\(J = 18\)) as the more complex one.
Figure 10: Scree plots of inter-item correlation matrix eigenvalues. A sharply dominant first eigenvalue indicates unidimensionality; multiple eigenvalues above 1 indicate multidimensionality.
What to notice. A clearly unidimensional instrument shows one eigenvalue well above the rest, with a sharp elbow. A multidimensional instrument shows multiple eigenvalues above 1, with a more gradual decline. The threshold is conventional rather than principled, which is what makes Holland’s result attractive in theory: rank is a sharper notion than “above 1”.
In practice the attraction does not survive contact with the estimator. The rank claim is about \(\boldsymbol{\Gamma}\)’s off-diagonals with the diagonal free, and Section 3 shows what it takes to test even \(D = 1\) — a weighted rank-1 fit to the off-diagonals plus a simulated reference distribution, because the eigenvalue count has no power and the share has no interpretable scale without one. Doing that for general \(D\) means minimum-rank fitting at each candidate rank with a reference distribution for each, on a model that needs all \(2^J\) cells. For the \(J\) at which instruments are actually screened for dimensionality, the SOE cannot be fitted at all. The correlation-matrix scree above is the pragmatic tool; the \(\boldsymbol{\Gamma}\) result is a statement about what dimensionality means rather than a procedure for measuring it.
Section 6: How many parameters does a long test support?
Holland (1990) establishes a conservation law: a \(D\)-dimensional test of \(J\) items supports at most \((D+1)J\) parameters as \(J \to \infty\). For \(D = 1\) this means at most \(2J\) parameters — exactly what the 2PL provides. The 3PL has \(3J\) parameters, exceeding the limit. Holland’s prediction is therefore that the marginal gain from 2PL → 3PL should shrink as \(J\) grows, and may become negative (overfitting) for long tests.
We test this using the InterModel Vigorish [IMV; see also Holland (1990)] — a log-score measure of predictive improvement — across IRW datasets spanning a range of test lengths.
Figure 11: IMV gain from 1PL → 2PL (blue) and 2PL → 3PL (red) as a function of test length J. Holland’s conservation law predicts the 2PL → 3PL gain shrinks toward zero as J grows.
What to notice. The 1PL → 2PL gain (blue) is typically positive and stable across \(J\): freeing discrimination improves fit regardless of test length, consistent with real variation in item discriminations. The 2PL → 3PL gain (red) should, if Holland is right, trend toward zero as \(J\) grows. If it remains positive at large \(J\), that suggests genuine multidimensionality or pervasive guessing. The pattern across the IRW is the empirical verdict.
WarningLimitations
Pattern sparsity. Sections 1–3 require \(J \leq 10\) for pattern enumeration to be feasible. The SOE loglinear approach does not generalise to large \(J\).
Normal ability assumption. The quadratic \(\log \pi_k\) result (Section 4) holds asymptotically under normality; skewed or multimodal ability distributions can mimic the pattern for the wrong reasons.
3PL instability. The 3PL is difficult to estimate for short tests and small samples. IMV values for 2PL → 3PL should be interpreted cautiously for \(J < 10\).
Dataset selection.irw_filter() conditions are heuristic. The “unidimensional” and “complex” candidates in Section 5 were selected on eigenvalue structure alone. Users with substantive knowledge of specific IRW tables should override the automatic selection.
These results were computed against approximately IRW v607 (the corpus as of October 22, 2026).
NoteHow this version was determined
The version was inferred from when the results were computed, not recorded during the run, so it may be off by a version or more. To pin data to an exact version, use irw_use_version().
Results were computed on October 22, 2026 with irw version 1.1.2. To regenerate:
# 1. Re-run the compute script (from the project root)source("vignettes/dutch_identity_compute.R")# 2. Re-render this pagequarto::quarto_render("vignettes/dutch_identity.qmd")
The model-fitting and Γ-diagnostic helpers live in vignettes/dutch_identity_helpers.R, which dutch_identity_compute.R sources. vignettes/dutch_identity_gamma_rerun.R updates only the Section 3 Γ diagnostics in the cached result, refitting the loglinear models from the stored pattern counts and asserting that Table 1 is unchanged, so a figure fix does not require re-running the corpus-wide sections.
Holland, Paul W. 1990. “The Dutch Identity: A New Tool for the Study of Item Response Models.”Psychometrika 55 (1): 5–18. https://doi.org/10.1007/BF02294739.
