Multidimensional IRT & the Factor-Analysis Bridge

A Discrimination Vector Is a Loading Vector

A 2PL measures one ability. Real instruments measure several correlated traits at once, so give each person a vector of latent traits and each item a vector of discriminations. At which point the model stops being only psychometrics: those discrimination vectors are factor loadings, and multidimensional IRT is confirmatory factor analysis for categorical items. The 2PL is the one-factor special case, and a discrimination has always been a loading — the project's title bridge, demonstrated rather than asserted.

Augmentation plus a parameter-expanded factor step

The sampler combines the Albert–Chib augmentation with the parameter-expanded factor Gibbs of Ghosh & Dunson: sample factors and loadings with an unconstrained factor covariance, which is conjugate inverse-Wishart, then standardise to unit-variance factors afterwards. That recovers the identified loadings and the factor correlation together, without a rotation step. A confirmatory mask fixes which items load on which trait, which is what identifies the solution in the first place.

Extraversion and neuroticism

On two Big Five scales — 2,617 respondents, five extraversion and five neuroticism items, binary-scored — all ten items load positively on their own trait and the two factors correlate −0.29, 95% CrI [−0.35, −0.24]: more extraverted respondents tend to be less neurotic. That is the familiar modest E–N link, recovered here directly from binary item responses. R's mirt gets −0.30 and PyMC −0.33, so the substantive finding is stable across three independent fits.

Checking the bridge

The bridge itself is checked explicitly. Refitting the five extraversion items alone as a one-factor model — which is precisely a 2PL — reproduces the two-factor loadings at correlation 0.970. Simple structure puts each item on one axis, which is what confirmatory factor analysis calls a clean factor pattern.

Two Loading Tables, One Fit

The same species of problem as the 2PL project's metric mismatch appears here, with a different transformation. The two engines print loading tables that look unrelated — 0.82, 1.33, 0.74, … from the sampler against 0.63, 0.80, 0.59, … from mirt — which is easy to take as two different results. They are the same numbers under the factor-analysis standardisation $\lambda=a/\sqrt{1+a^2}$, which converts a normal-ogive discrimination into a loading on a unit-variance latent response. Applied across all ten items the maximum absolute difference is 0.011 and the correlation 0.9991 — agreement to rounding. Both notebooks perform the conversion rather than comparing orderings by eye.

Why two packages disagree about a factor correlation

A second comparison deserves the same care. R's psych::fa reports a factor correlation of −0.22 against mirt's −0.30 — a third apart, and the loading table explains why: psych is fitted exploratory with an oblimin rotation, so it is free to place cross-loadings — and it does, N4 picking up −0.35 on the extraversion factor, N5 −0.14, E5 +0.13. Every bit of between-factor association those cross-paths absorb is association the factor correlation no longer has to carry. mirt is confirmatory, with those paths fixed at zero, so the same covariance has one place left to go. The qualitative conclusion is robust; the magnitude of a factor correlation depends on what else the model was allowed to explain, which is a general caution about reading one off any rotated solution.

Notebooks

Downloads

Item-Factor Module — Source Code

"""
factorirt.py -- MULTIDIMENSIONAL IRT = the item factor-analysis model (from scratch).

Backs the notebooks in  "Multidimensional IRT & the Factor-Analysis Bridge".

A 2PL measures one ability. Real instruments often measure SEVERAL correlated traits at once
(extraversion AND neuroticism, verbal AND quantitative). Multidimensional IRT gives each person
a VECTOR of latent traits theta_i and each item a vector of DISCRIMINATIONS -- which are exactly
FACTOR LOADINGS. In the normal-ogive form,

    P(x_ij = 1 | theta_i) = Phi( a_j' theta_i + d_j ),     theta_i ~ N(0, Sigma),

a_j is item j's loading vector, d_j the easiness, and Sigma the correlation between the latent
factors. This IS confirmatory factor analysis for binary items (item factor analysis): the 2PL
is the one-factor special case, and the discrimination is a loading. A CONFIRMATORY model fixes
which items load on which factor (simple structure), which identifies the loadings without
rotation.

From scratch it is the ALBERT-CHIB augmentation again -- draw z_ij ~ N(a_j'theta_i + d_j, 1)
truncated by the sign of x_ij -- combined with the PARAMETER-EXPANDED Gibbs for factor analysis
(Ghosh & Dunson 2009): sample the factors and loadings with an UNCONSTRAINED factor covariance
(inverse-Wishart, conjugate), then STANDARDISE afterwards to unit-variance factors, which recovers
the identified loadings and the factor correlation. Per-factor sign is fixed by keeping each
factor's loadings summing positive.
"""

import numpy as np
from scipy.special import ndtr as Phi, ndtri as Phinv
from scipy.stats import invwishart


def _rtrunc_sign(mean, pos, rng):
    lo = np.where(pos, Phi(-mean), 0.0); hi = np.where(pos, 1.0, Phi(-mean))
    u = lo + rng.random(mean.shape) * (hi - lo)
    return mean + Phinv(np.clip(u, 1e-12, 1 - 1e-12))


# --------------------------------------------------------------------------- #
#  simulation                                                                   #
# --------------------------------------------------------------------------- #

def simulate_mirt(N, A, d, rng, Sigma=None):
    """A:(J,D) loadings, d:(J,) intercepts, Sigma:(D,D) factor correlation."""
    A = np.asarray(A, float); d = np.asarray(d, float); J, D = A.shape
    if Sigma is None:
        Sigma = np.eye(D)
    theta = rng.multivariate_normal(np.zeros(D), Sigma, N)
    z = theta @ A.T + d + rng.standard_normal((N, J))
    return (z > 0).astype(float), theta


# --------------------------------------------------------------------------- #
#  Gibbs sampler (Albert-Chib augmentation + parameter-expanded factor Gibbs)   #
# --------------------------------------------------------------------------- #

def mirt_gibbs(X, mask, rng, draws=2500, burn=1500, load_sd=2.0):
    """multidimensional 2PL / item factor analysis. mask:(J,D) boolean -- which factors each item
    loads on (a confirmatory simple structure). Returns standardised loadings A:(draws,J,D),
    intercepts d:(draws,J), factor correlations R:(draws,D,D) and factors theta:(draws,N,D)."""
    X = np.asarray(X, float); N, J = X.shape; D = mask.shape[1]
    A = mask.astype(float) * 0.7; d = np.zeros(J); theta = rng.standard_normal((N, D)); Sig = np.eye(D)
    AS = np.empty((draws, J, D)); DD = np.empty((draws, J)); RR = np.empty((draws, D, D)); TH = np.empty((draws, N, D))
    for it in range(draws + burn):
        z = _rtrunc_sign(theta @ A.T + d, X > 0.5, rng)
        # theta_i | z, A, d, Sig
        Vinv = A.T @ A + np.linalg.inv(Sig); V = np.linalg.inv(Vinv); Lv = np.linalg.cholesky(V)
        M = ((z - d) @ A) @ V.T
        theta = M + rng.standard_normal((N, D)) @ Lv.T
        # (a_j, d_j) | z, theta   -- masked regression of z_j on [theta_cols, 1]
        for j in range(J):
            cols = np.where(mask[j])[0]; p = len(cols)
            Xd = np.column_stack([theta[:, cols], np.ones(N)])
            prior = np.diag([1.0 / load_sd ** 2] * p + [1e-6])
            Vj = np.linalg.inv(Xd.T @ Xd + prior); mj = Vj @ (Xd.T @ z[:, j])
            dr = mj + np.linalg.cholesky(Vj) @ rng.standard_normal(p + 1)
            A[j, cols] = dr[:p]; d[j] = dr[-1]
        for k in range(D):                                # per-factor sign identification
            if A[:, k].sum() < 0:
                A[:, k] = -A[:, k]; theta[:, k] = -theta[:, k]
        # Sigma | theta  (unconstrained inverse-Wishart, then standardise to a correlation)
        Sig = np.atleast_2d(invwishart.rvs(df=N + D + 1, scale=np.eye(D) + theta.T @ theta, random_state=rng))
        s = np.sqrt(np.diag(Sig))
        if it >= burn:
            i = it - burn
            AS[i] = A * s[None, :]; DD[i] = d; RR[i] = Sig / np.outer(s, s); TH[i] = theta / s[None, :]
    return dict(A=AS, d=DD, R=RR, theta=TH)

References

• Bock, R. D., Gibbons, R. & Muraki, E. (1988). Full-information item factor analysis. Applied Psychological Measurement 12(3), 261–280. — the model this project fits, and the name of the bridge
• Béguin, A. A. & Glas, C. A. W. (2001). MCMC estimation and some model-fit analysis of multidimensional IRT models. Psychometrika 66(4), 541–561. — the augmentation sampler in the multidimensional case
• Ghosh, J. & Dunson, D. B. (2009). Default prior distributions and efficient posterior computation in Bayesian factor analysis. JCGS 18(2), 306–320. — the parameter expansion that makes the factor covariance conjugate
• Takane, Y. & de Leeuw, J. (1987). On the relationship between item response theory and factor analysis of discretized variables. Psychometrika 52(3), 393–408. — the formal equivalence, including the standardisation used here
• Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. — the exploratory tetrachoric factor analysis compared against the confirmatory fit
• Chalmers, R. P. (2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software 48(6). — the confirmatory reference fit