Multidimensional IRT & the Factor-Analysis Bridge¶
Discriminations are factor loadings¶
A 2PL measures one ability. Real instruments usually tap several correlated traits — a Big-Five questionnaire measures extraversion and neuroticism, an exam measures verbal and quantitative skill. Multidimensional IRT gives each person a vector of latent traits $\theta_i$ and each item a vector of discriminations, which are exactly factor loadings: $$P(x_{ij}=1\mid\theta_i)=\Phi(a_j^\top\theta_i+d_j),\qquad \theta_i\sim N(0,\Sigma).$$ Here $a_j$ is item $j$'s loading vector, $d_j$ its easiness, and $\Sigma$ the correlation between the latent factors. This is confirmatory factor analysis for binary items (item factor analysis) — and the 2PL is the one-factor special case, where the single discrimination is a single loading. A confirmatory model fixes which items load on which factor ("simple structure"), which identifies the loadings without any rotation.
From scratch it is the Albert–Chib augmentation once more, combined with the parameter-expanded factor-analysis Gibbs (Ghosh–Dunson 2009): sample factors and loadings with an unconstrained factor covariance (inverse-Wishart, conjugate), then standardise to unit-variance factors to recover the identified loadings and the factor correlation. We validate recovery, fit two Big-Five scales — extraversion and neuroticism — as a two-factor model, read the loadings and the inter-factor correlation, connect it to the 2PL and to classical factor analysis, and cross-check in PyMC.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import factorirt as F
rng = np.random.default_rng(3)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("Multidimensional IRT = item factor analysis: a discrimination vector is a loading vector.")
Multidimensional IRT = item factor analysis: a discrimination vector is a loading vector.
1. Does it work? — recovering loadings and the factor correlation¶
Simulated responses of 3000 people to ten binary items: five load on factor 1, five on factor 2 (a confirmatory simple structure), and the two factors are correlated ($\rho=0.4$). The augmentation + parameter-expanded Gibbs should recover the loadings and the correlation.
J,D=10,2; A_t=np.zeros((J,D)); A_t[:5,0]=[1.2,1.0,0.9,1.4,0.8]; A_t[5:,1]=[1.1,0.9,1.3,0.7,1.0]
rho=0.4; Sig=np.array([[1,rho],[rho,1]]); X,th=F.simulate_mirt(3000,A_t,np.zeros(J),rng,Sig)
mask=np.zeros((J,D),bool); mask[:5,0]=True; mask[5:,1]=True
r=F.mirt_gibbs(X,mask,rng,draws=1200,burn=800); A_h=r["A"].mean(0); R_h=r["R"].mean(0)
lt=np.r_[A_t[:5,0],A_t[5:,1]]; lh=np.r_[A_h[:5,0],A_h[5:,1]]
fig,ax=plt.subplots(1,2,figsize=(11,4.2))
ax[0].scatter(lt,lh,color=BLUE); ax[0].plot([0,1.6],[0,1.6],"k--",lw=1); ax[0].set_xlabel("true loading"); ax[0].set_ylabel("estimated loading"); ax[0].set_title(f"loadings (r={np.corrcoef(lt,lh)[0,1]:.2f}, RMSE={np.sqrt(np.mean((lt-lh)**2)):.2f})")
rd=r["R"][:,0,1]; ax[1].hist(rd,bins=30,color=PURP,alpha=.7); ax[1].axvline(rho,color="k",ls="--",lw=1.5,label=f"true ρ={rho}"); ax[1].axvline(rd.mean(),color=RED,lw=1.5,label=f"estimated {rd.mean():.2f}")
ax[1].set_xlabel("factor correlation"); ax[1].set_ylabel("posterior"); ax[1].set_title("Inter-factor correlation"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("Loadings and the factor correlation are recovered -- the confirmatory simple structure identifies them without")
print("rotation, and the parameter-expanded step gives the correlation between the two latent traits.")
Loadings and the factor correlation are recovered -- the confirmatory simple structure identifies them without rotation, and the parameter-expanded step gives the correlation between the two latent traits.
2. Two Big-Five scales — extraversion and neuroticism¶
2617 respondents answered five extraversion items (E1–E5) and five neuroticism items (N1–N5), scored agree/disagree. We fit a two-factor model with the natural simple structure — E items load on an extraversion factor, N items on a neuroticism factor — and let the two factors correlate. The loadings say how strongly each item marks its trait; the factor correlation is a substantive personality finding.
d=pd.read_csv("bfi_EN.csv"); X=d.to_numpy().astype(float); items=list(d.columns)
mask=np.zeros((10,2),bool); mask[:5,0]=True; mask[5:,1]=True
rb=F.mirt_gibbs(X,mask,rng,draws=2500,burn=1500); A=rb["A"].mean(0); R=rb["R"].mean(0); rho_lo,rho_hi=np.percentile(rb["R"][:,0,1],[2.5,97.5])
load=np.where(mask, A, np.nan)
print("factor loadings (blank = fixed to zero by the confirmatory structure):")
print(pd.DataFrame({"Extraversion":load[:,0].round(2),"Neuroticism":load[:,1].round(2)}, index=items).to_string(na_rep=""))
print(f"\nExtraversion-Neuroticism factor correlation = {R[0,1]:.2f} 95% CrI [{rho_lo:.2f}, {rho_hi:.2f}]")
fig,ax=plt.subplots(1,2,figsize=(12,4.3))
y=np.arange(10); ax[0].barh(y[:5], A[:5,0], color=BLUE, label="Extraversion"); ax[0].barh(y[5:], A[5:,1], color=RED, label="Neuroticism")
ax[0].set_yticks(y); ax[0].set_yticklabels(items); ax[0].invert_yaxis(); ax[0].set_xlabel("loading on own factor"); ax[0].set_title("Item loadings"); ax[0].legend(frameon=False)
tp=rb["theta"].mean(0); ax[1].scatter(tp[:,0],tp[:,1],s=5,color=GREY,alpha=.25)
ax[1].set_xlabel("Extraversion factor"); ax[1].set_ylabel("Neuroticism factor"); ax[1].set_title(f"Estimated factor scores (corr {np.corrcoef(tp[:,0],tp[:,1])[0,1]:.2f})")
plt.tight_layout(); plt.show()
a_flat = A.sum(1) # simple structure: one nonzero loading per item
lam = a_flat / np.sqrt(1 + a_flat**2) # normal-ogive discrimination -> standardised loading
lam_mirt = np.array([0.63,0.80,0.59,0.77,0.60,0.84,0.84,0.78,0.62,0.57]) # mirt, standardised
print("\nWHICH SCALE. These are DISCRIMINATIONS on the normal-ogive scale; mirt prints STANDARDISED")
print("loadings, i.e. correlations between an item's latent response and its factor. The two are one")
print("quantity under lambda = a / sqrt(1 + a^2) -- the factor-analysis normalisation that divides by")
print("the latent response's total standard deviation:\n")
print(" item a a/sqrt(1+a^2) mirt")
for nm, ai, li, mi in zip(items, a_flat, lam, lam_mirt):
print(" %-6s %6.2f %13.2f %9.2f" % (nm, ai, li, mi))
print("\n max absolute difference %.3f, correlation %.4f -- agreement to rounding on all ten items."
% (np.abs(lam - lam_mirt).max(), np.corrcoef(lam, lam_mirt)[0,1]))
print("Two tables that look unrelated at a glance are the same fit in two conventions. Worth converting")
print("explicitly rather than eyeballing the orderings and calling it agreement.\n")
print(f"All ten items load positively on their own trait. The two factors correlate NEGATIVELY (rho={R[0,1]:.2f}): more")
print("extraverted respondents tend to be less neurotic -- the well-known modest E-N link in the Big Five, recovered")
print("directly from binary item responses by an item factor analysis.")
factor loadings (blank = fixed to zero by the confirmatory structure):
Extraversion Neuroticism
E1 0.82
E2 1.33
E3 0.74
E4 1.18
E5 0.74
N1 1.53
N2 1.54
N3 1.27
N4 0.80
N5 0.71
Extraversion-Neuroticism factor correlation = -0.29 95% CrI [-0.35, -0.24]
WHICH SCALE. These are DISCRIMINATIONS on the normal-ogive scale; mirt prints STANDARDISED loadings, i.e. correlations between an item's latent response and its factor. The two are one quantity under lambda = a / sqrt(1 + a^2) -- the factor-analysis normalisation that divides by the latent response's total standard deviation: item a a/sqrt(1+a^2) mirt E1 0.82 0.63 0.63 E2 1.33 0.80 0.80 E3 0.74 0.59 0.59 E4 1.18 0.76 0.77 E5 0.74 0.60 0.60 N1 1.53 0.84 0.84 N2 1.54 0.84 0.84 N3 1.27 0.79 0.78 N4 0.80 0.63 0.62 N5 0.71 0.58 0.57 max absolute difference 0.011, correlation 0.9991 -- agreement to rounding on all ten items. Two tables that look unrelated at a glance are the same fit in two conventions. Worth converting explicitly rather than eyeballing the orderings and calling it agreement. All ten items load positively on their own trait. The two factors correlate NEGATIVELY (rho=-0.29): more extraverted respondents tend to be less neurotic -- the well-known modest E-N link in the Big Five, recovered directly from binary item responses by an item factor analysis.
Reading a loading, a communality and a factor correlation¶
The table above prints discriminations; mirt prints standardised loadings. They
are the same fit, and $\lambda_j = a_j/\sqrt{1+a_j^2}$ converts between them — but only the
standardised form has a direct reading, because it is a correlation: the correlation between the
item's latent response and the factor it loads on.
That makes its square immediately meaningful. $\lambda_j^2$ is the item's communality — the share of its latent-response variance the trait accounts for, with the remainder being item-specific. N1 and N2 at $\lambda = 0.84$ have communalities of about 71%, so most of what moves those responses is neuroticism itself. N5 at $\lambda = 0.58$ sits at about 34%: still a genuine indicator, but two-thirds of its variation is something other than the trait. That gap is the practical difference between a good item and a passable one, and it is invisible in the raw discriminations, where 1.53 against 0.71 looks like a factor of two rather than a doubling of explained variance.
The factor correlation of $-0.29$ is a correlation between the latent traits, not between observed scale scores. This matters: an observed correlation is attenuated by measurement error and will be smaller, so the two numbers are not comparable and the latent one is the quantity a substantive claim should be about. Its magnitude is modest — squared, extraversion and neuroticism share about 8% of their variance — and the interval $[-0.35, -0.24]$ excludes zero, so the sign is established while the effect stays small. "Modest negative link" is the honest phrasing; "the two traits are opposites" is not.
Two structural points about what the model was and was not free to find. The blanks in the loading table are fixed at zero by the confirmatory design, not estimated and found to be small — the model was told each item belongs to one factor and asked how strongly, so a clean pattern here is not evidence that no cross-loadings exist. And because that simple structure is imposed, the loadings are identified without rotation: an unconstrained factor model can be spun to infinitely many equally-fitting solutions, which is why exploratory analyses must choose a rotation and confirmatory ones need not. The parameter-expanded inverse-Wishart step is what lets the sampler draw a correlation matrix — unit diagonal, which is an awkward constraint to sample directly — by working with an unconstrained covariance and rescaling it, avoiding the slow mixing a directly constrained sampler would suffer.
3. The bridge — the 2PL is a one-factor model, and this is factor analysis¶
Two connections make the picture whole. First, fitting one factor to the extraversion items alone returns exactly a 2PL: the loadings are the discriminations. Second, the whole model is factor analysis of binary items — the same latent-variable idea as classical factor analysis, but through a probit link, so the loadings are on the familiar factor-analytic scale.
XE=X[:,:5]; mask1=np.ones((5,1),bool)
r1=F.mirt_gibbs(XE,mask1,rng,draws=2000,burn=1200); a1=r1["A"].mean(0)[:,0]
print("Extraversion items, one-factor model (= a 2PL):")
print(pd.DataFrame({"1-factor loading":a1.round(2),"loading in 2-factor model":A[:5,0].round(2)}, index=items[:5]).to_string())
print(f"\nagreement: {np.corrcoef(a1,A[:5,0])[0,1]:.3f} -- the one-factor loadings are the item discriminations; MIRT is")
print("multidimensional 2PL, and a single dimension is the ordinary 2PL. The discrimination has always been a loading.")
# a picture: two-factor loading space
fig,ax=plt.subplots(figsize=(6.4,5))
ax.scatter(A[:5,0], np.zeros(5), color=BLUE, s=60); ax.scatter(np.zeros(5), A[5:,1], color=RED, s=60)
for i in range(5): ax.annotate(items[i],(A[i,0],0),fontsize=8,xytext=(0,6),textcoords="offset points")
for i in range(5): ax.annotate(items[5+i],(0,A[5+i,1]),fontsize=8,xytext=(6,0),textcoords="offset points")
ax.axhline(0,color=GREY,lw=.8); ax.axvline(0,color=GREY,lw=.8)
ax.set_xlabel("loading on Extraversion"); ax.set_ylabel("loading on Neuroticism"); ax.set_title("Confirmatory simple structure: each item on one axis")
plt.tight_layout(); plt.show()
print("Simple structure puts each item on one axis -- what confirmatory factor analysis calls a clean factor pattern.")
Extraversion items, one-factor model (= a 2PL):
1-factor loading loading in 2-factor model
E1 0.85 0.82
E2 1.17 1.33
E3 0.78 0.74
E4 1.20 1.18
E5 0.76 0.74
agreement: 0.970 -- the one-factor loadings are the item discriminations; MIRT is
multidimensional 2PL, and a single dimension is the ordinary 2PL. The discrimination has always been a loading.
Simple structure puts each item on one axis -- what confirmatory factor analysis calls a clean factor pattern.
4. Cross-check in PyMC¶
The multidimensional 2PL sampled directly: standard-normal factors with a free correlation $\rho$, positive loadings on each item's own factor, and a Bernoulli likelihood through the probit link. NUTS samples the factor scores and the loadings; we compare loadings and the E–N correlation with the from-scratch sampler. (We fit a representative subsample for speed.)
import pymc as pm, pytensor.tensor as pt
sub=rng.choice(len(X),1200,replace=False); Xs=X[sub]; Ns=len(Xs)
with pm.Model() as mod:
rho=pm.Uniform("rho",-0.95,0.95)
one=pt.ones(())
cov=pt.stack([pt.stack([one,rho]), pt.stack([rho,one])])
theta=pm.MvNormal("theta", mu=np.zeros(2), cov=cov, shape=(Ns,2))
aE=pm.HalfNormal("aE",2.0,shape=5); aN=pm.HalfNormal("aN",2.0,shape=5); d_=pm.Normal("d",0,3,shape=10)
Amat=pt.concatenate([pt.stack([aE,pt.zeros(5)],axis=1), pt.stack([pt.zeros(5),aN],axis=1)],axis=0) # (10,2)
eta=pt.dot(theta,Amat.T)+d_
pm.Bernoulli("x", p=pm.math.invprobit(eta), observed=Xs)
idata=pm.sample(500, tune=1000, chains=4, target_accept=0.9, random_seed=4, progressbar=False)
aE_pm=idata.posterior["aE"].mean(("chain","draw")).values; aN_pm=idata.posterior["aN"].mean(("chain","draw")).values
rho_pm=float(idata.posterior["rho"].mean())
print(f"factor correlation: from-scratch {R[0,1]:.2f} PyMC {rho_pm:.2f}")
print("loadings (from-scratch / PyMC):")
for i in range(5): print(f" {items[i]:3s} E {A[i,0]:.2f} / {aE_pm[i]:.2f} {items[5+i]:3s} N {A[5+i,1]:.2f} / {aN_pm[i]:.2f}")
fig,ax=plt.subplots(figsize=(6.6,4)); ax.scatter(np.r_[A[:5,0],A[5:,1]], np.r_[aE_pm,aN_pm], color=BLUE)
ax.plot([0,1.7],[0,1.7],"k--",lw=1); ax.set_xlabel("from-scratch loading"); ax.set_ylabel("PyMC loading"); ax.set_title("Loadings: from-scratch vs PyMC")
plt.tight_layout(); plt.show()
print("Same loadings, same negative factor correlation: the augmentation Gibbs and PyMC's direct NUTS agree on the")
print("two-factor structure and on the extraversion-neuroticism link.")
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [rho, theta, aE, aN, d]
Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 24 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
factor correlation: from-scratch -0.29 PyMC -0.33 loadings (from-scratch / PyMC): E1 E 0.82 / 0.88 N1 N 1.53 / 1.55 E2 E 1.33 / 1.26 N2 N 1.54 / 1.63 E3 E 0.74 / 0.64 N3 N 1.27 / 1.25 E4 E 1.18 / 1.10 N4 N 0.80 / 0.75 E5 E 0.74 / 0.71 N5 N 0.71 / 0.67
Same loadings, same negative factor correlation: the augmentation Gibbs and PyMC's direct NUTS agree on the two-factor structure and on the extraversion-neuroticism link.
5. Summary¶
Multidimensional IRT is item factor analysis: each item's discrimination becomes a loading vector, and the model estimates how the items load on several correlated latent traits. From scratch it is the Albert–Chib augmentation plus a parameter-expanded factor-analysis Gibbs — recovering loadings and the factor correlation from simulated data, and, on the Big Five extraversion and neuroticism scales, finding that all items load cleanly on their own factor and that the two traits are modestly negatively correlated ($\rho\approx-0.3$) — the classic personality result, obtained directly from binary responses and confirmed in PyMC.
This fills the factor-analysis gap in the portfolio and closes several loops: the 2PL is the one-factor case (a discrimination is a loading); the probit-factor engine is the same one the Conditional-Dependence LCA used to put a continuous trait inside latent classes; and the augmentation is Albert (1992), shared with the Tobit, multivariate-probit and probit-LCA notebooks. Confirmatory simple structure is exactly what CFA imposes; the frequentist item-factor-analysis fits (mirt, psych) are in the companion R notebook. Next in the arc: differential item functioning and explanatory IRT.