Multidimensional IRT & the Factor-Analysis Bridge (R)¶
mirt (confirmatory) and psych (item factor analysis)¶
The R counterpart to mirt_python.ipynb. Where the Python notebook builds the multidimensional 2PL / item factor analysis from scratch (Albert–Chib augmentation + parameter-expanded factor Gibbs) and cross-checks in PyMC, this notebook uses the field-standard packages: mirt fits a confirmatory two-factor IRT model directly, and psych performs the classical factor analysis of tetrachoric correlations — the two views of the same item factor analysis. We fit the Big-Five extraversion and neuroticism scales (five binary items each) and read the loadings and the inter-factor correlation.
options(repr.plot.width=12, repr.plot.height=4.6)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages({library(mirt); library(psych)})
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; GREY<-"#718096"
X<-read.csv("bfi_EN.csv"); cat(nrow(X), "respondents x", ncol(X), "binary items (E1-E5, N1-N5)\n")
Warning message: "package 'mirt' was built under R version 4.6.1"
Warning message: "package 'psych' was built under R version 4.6.1"
2617 respondents x 10 binary items (E1-E5, N1-N5)
1. Confirmatory two-factor IRT with mirt¶
mirt.model specifies the simple structure — E items load on an Extraversion factor, N items on a Neuroticism factor, and the two factors covary. mirt estimates the loadings (standardised) and the factor correlation by marginal maximum likelihood.
mod <- mirt.model('Extra = 1-5
Neuro = 6-10
COV = Extra*Neuro')
m <- mirt(X, mod, itemtype="2PL", verbose=FALSE)
sm <- summary(m, verbose=FALSE)
cat("standardised factor loadings:\n"); print(round(sm$rotF, 2))
fc <- sm$fcor
cat(sprintf("\nExtraversion-Neuroticism factor correlation = %.2f\n", fc[1,2]))
L <- sm$rotF
par(mar=c(4,4,3,1)); y<-1:10
barplot(t(cbind(ifelse(L[,1]>L[,2],L[,1],0), ifelse(L[,2]>=L[,1],L[,2],0))), beside=FALSE, horiz=TRUE,
col=c(BLUE,RED), names.arg=rownames(L), las=1, xlab="loading on own factor", main="mirt confirmatory loadings")
legend("bottomright", c("Extraversion","Neuroticism"), fill=c(BLUE,RED), bty="n")
cat("\nWHICH SCALE. mirt prints STANDARDISED loadings -- correlations between an item's latent\n")
cat("response and its factor. The from-scratch sampler in the Python notebook prints DISCRIMINATIONS on\n")
cat("the normal-ogive scale (0.82, 1.33, ...), which look unrelated until converted by lambda = a/sqrt(1+a^2).\n")
cat("Applying that to all ten items reproduces the column above to within 0.011, so the two engines agree\n")
cat("exactly -- the appearance of disagreement is entirely a difference of convention.\n\n")
cat("All ten items load cleanly on their own factor; the two factors correlate negatively -- more extraverted\n")
cat("respondents tend to be less neurotic, the classic Big-Five link, and the same result as the from-scratch fit.\n")
standardised factor loadings:
Extra Neuro E1 0.63 0.00 E2 0.80 0.00 E3 0.59 0.00 E4 0.77 0.00 E5 0.60 0.00 N1 0.00 0.84 N2 0.00 0.84 N3 0.00 0.78 N4 0.00 0.62 N5 0.00 0.57
Extraversion-Neuroticism factor correlation = -0.30
WHICH SCALE. mirt prints STANDARDISED loadings -- correlations between an item's latent
response and its factor. The from-scratch sampler in the Python notebook prints DISCRIMINATIONS on
the normal-ogive scale (0.82, 1.33, ...), which look unrelated until converted by lambda = a/sqrt(1+a^2).
Applying that to all ten items reproduces the column above to within 0.011, so the two engines agree
exactly -- the appearance of disagreement is entirely a difference of convention.
All ten items load cleanly on their own factor; the two factors correlate negatively -- more extraverted
respondents tend to be less neurotic, the classic Big-Five link, and the same result as the from-scratch fit.
2. Classical item factor analysis with psych¶
The same model seen the factor-analytic way: compute the tetrachoric correlations of the binary items (the correlations of the underlying continuous propensities) and factor-analyse them with an oblique (correlated-factors) rotation. psych::fa should give the same loading pattern and the same factor correlation as mirt — because multidimensional IRT is factor analysis of the tetrachoric correlations.
fa2 <- fa(X, nfactors=2, cor="tet", rotate="oblimin", fm="ml")
cat("psych fa loadings (oblimin):\n"); print(round(unclass(fa2$loadings), 2))
cat(sprintf("\npsych factor correlation = %.2f (mirt gave %.2f)\n", fa2$Phi[1,2], fc[1,2]))
# align the two factors and compare loadings
mir_load <- apply(sm$rotF, 1, max) # each item's loading on its own factor (mirt)
fa_load <- apply(abs(unclass(fa2$loadings)), 1, max)
par(mar=c(4,4,3,1)); plot(mir_load, fa_load, pch=19, col=BLUE, xlab="mirt loading", ylab="psych fa loading", main=sprintf("Loadings agree (r = %.3f)", cor(mir_load, fa_load)))
text(mir_load, fa_load, rownames(sm$rotF), pos=3, cex=.7); abline(0,1,lty=2,col=GREY)
cat("\nThe two factor correlations are NOT the same number: psych gives -0.22 against mirt's -0.30, a gap\n")
cat("of about a third. That is a difference of specification rather than of data, and the loading table above\n")
cat("shows why: psych is fitted EXPLORATORY with an oblimin rotation, so it is free to place cross-loadings,\n")
cat("and it does -- N4 picks up -0.35 on the extraversion factor, N5 -0.14, E5 +0.13. Every bit of\n")
cat("between-factor association those cross-paths absorb is association the factor correlation no longer has\n")
cat("to carry. mirt is CONFIRMATORY: those paths are fixed at zero, so the covariance has one place to go.\n")
cat("\nThe robust statement is the qualitative one -- multidimensional IRT and item factor analysis are the\n")
cat("same model, both recover the classic negative extraversion-neuroticism link, and the loadings track each\n")
cat("other closely. The MAGNITUDE of a factor correlation depends on what else the model was allowed to\n")
cat("explain -- a general caution about reading one off any rotated solution.\n")
psych fa loadings (oblimin):
ML1 ML2 E1 0.09 0.67 E2 -0.13 0.75 E3 0.06 0.60 E4 -0.03 0.73 E5 0.13 0.65 N1 0.90 0.07 N2 0.88 0.02 N3 0.73 -0.05 N4 0.48 -0.35 N5 0.50 -0.14
psych factor correlation = -0.22 (mirt gave -0.30)
The two factor correlations are NOT the same number: psych gives -0.22 against mirt's -0.30, a gap
of about a third. That is a difference of specification rather than of data, and the loading table above
shows why: psych is fitted EXPLORATORY with an oblimin rotation, so it is free to place cross-loadings,
and it does -- N4 picks 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: those paths are fixed at zero, so the covariance has one place to go.
The robust statement is the qualitative one -- multidimensional IRT and item factor analysis are the
same model, both recover the classic negative extraversion-neuroticism link, and the loadings track each
other closely. The MAGNITUDE of a factor correlation depends on what else the model was allowed to
explain -- a general caution about reading one off any rotated solution.
3. Summary¶
The standard R packages reproduce the Python result by two routes to the same item factor analysis: mirt fits the confirmatory two-factor IRT model directly (loadings + factor correlation by marginal ML), and psych::fa factor-analyses the tetrachoric correlations — and they agree on the loadings and on the negative extraversion–neuroticism correlation, matching the from-scratch Albert–Chib / parameter-expanded Gibbs and the PyMC fit in mirt_python.ipynb.
This is the factor-analysis bridge: a discrimination is a loading, the 2PL is the one-factor case, and confirmatory simple structure is exactly what confirmatory factor analysis imposes. mirt (confirmatory MIRT) and psych/lavaan (item factor analysis, CFA) are the field-standard tools, here the frequentist cross-check on the Bayesian from-scratch model. Next in the arc: differential item functioning and explanatory IRT.