Factor analysis II: confirmatory models, omega, and IRT as factor analysis

Module: Factor analysis
Tranche: core (a first course)
Before this: Factor analysis I: exploring dimensionality · From the 1PL to the 4PL
Builds toward: Dimensionality and multidimensional IRT · Measurement invariance: comparing groups with multigroup CFA · Structural equation modeling with lavaan · Continuous responses: sliders and bounded scales
Starts threads: Model-based reliability (omega) in place of alpha (→ Dimensionality and multidimensional IRT); Ordinal factor analysis is the graded response model (→ Dimensionality and multidimensional IRT · Continuous responses: sliders and bounded scales); Choose the estimator from the data: categories, shape, missingness (→ Measurement invariance: comparing groups with multigroup CFA · Continuous responses: sliders and bounded scales)
Picks up threads: Alpha is a lower bound on reliability (from Classical test theory and reliability); One factor accounts for the item correlations (from Factor analysis I: exploring dimensionality); The data fix the correlations, not the frame of the factors (from Factor analysis I: exploring dimensionality)
IRW tables: grit, neurodegenerative_huizinga_2019_dass
From EDUC 252: n/a
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.
Feedback: This page improved thanks to feedback from Lijin Zhang. All remaining errors are the responsibility of the page’s creator (BD).
Comments: open an issue on GitHub about this lesson.

Code

import {grid, palette} from "./widgets/irt.js"
import {zeros, implied, fitCFA, fitIndices, omegas, pOgive, ogiveToIRT, p2plAB} from "./widgets/fa-confirmatory.js"
import {quiz, predict} from "./widgets/quiz.js"

What this is for

The Grit-O scale was written around two facets of grit, perseverance of effort and consistency of interest. Does it measure one thing or two? Factor analysis I let the data propose a structure. This lesson states a structure first and asks whether the data agree: confirmatory factor analysis (CFA). A written-down structure also answers three more questions. How well does it fit? How reliable is a score built on it (omega)? And what is it, seen from item response theory? Each fit also needs an estimator, and we choose one out loud for both of our scales. Along the way we settle whether the DASS-21 yields one score or three.

Goals

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

write a confirmatory factor model in lavaan syntax, say how it is identified, and count its degrees of freedom;
read χ², CFI, TLI, RMSEA and SRMR, say what each responds to, and compare nested models;
compute omega and omega hierarchical, and say when omega equals alpha;
convert an ordinal factor loading and threshold into a slope and a difficulty, and check the conversion against mirt;
say what each estimator (ML, MLR, WLSMV) and each treatment of missing responses assumes, and what in the data bears on those assumptions.

Core ideas

Stating the structure first

A confirmatory model is the common factor model with most loadings fixed at zero before the data are seen. For grit, perseverance items load on one factor and consistency items on the other, with each item’s other loading fixed at 0. Jöreskog (1969) showed how to fit such a model by maximum likelihood and test it. We use the lavaan package (Rosseel, 2012), whose syntax reads almost like a path diagram. =~ means “is measured by”, and an item left off a line has its loading on that factor fixed at zero:

two <- "per =~ g1 + g4 + g6 + g9 + g10 + g12
        con =~ g2 + g3 + g5 + g7 + g8 + g11"
fit <- cfa(two, data = grit)

The factors may correlate (lavaan’s default), and each item keeps a unique variance.

Recall from Factor analysis I: with two or more factors, any rotation reproduces the correlations exactly as well as any other. The data fix the correlations, not the frame of the factors.

The zeros pin the frame, since rotating would move them. One indeterminacy is left: a factor has no natural unit. The model is identified once each factor’s scale is set, by fixing one loading per factor to 1 (lavaan’s default) or each factor’s variance to 1 (std.lv = TRUE). The two choices fit identically and report loadings in different units, as the grit output shows.

What can the data test? They supply \(p(p+1)/2\) variances and covariances. The model spends one parameter on each free loading, unique variance and factor correlation. The degrees of freedom are the difference, and each one is a constraint the data could violate.

Code

quiz({
  q: "Twelve items, two factors, each item loading on one factor only, factor variances fixed at 1 and the factors correlated. How many degrees of freedom does the model have?",
  options: [
    {text: "53", correct: true, why: "There are $12 \\times 13/2 = 78$ variances and covariances. The model estimates 12 loadings, 12 unique variances and 1 factor correlation, 25 in all: $78 - 25 = 53$."},
    {text: "54", correct: false, why: "That is the one-factor model, with 24 parameters. The factor correlation costs one more."},
    {text: "66", correct: false, why: "That counts only the 66 covariances and forgets the 12 variances the data also supply."}
  ]
})

Fit indices, and what they can’t tell you

Every fit index compares the covariance matrix the model implies with the one the data supply. The likelihood-ratio χ² tests exact fit: it is \(N - 1\) times the discrepancy, so a large sample rejects any model that is slightly wrong, and every model is. The other indices ask how wrong. RMSEA is the misfit per degree of freedom (Steiger, 1990). CFI (Bentler, 1990) and TLI (Tucker & Lewis, 1973) compare the model’s misfit with that of a baseline in which the items are uncorrelated. SRMR is the root mean square residual correlation.

Hu and Bentler (1999) proposed, from simulations of maximum-likelihood estimation, cutoffs near 0.95 for CFI and TLI, 0.06 for RMSEA and 0.08 for SRMR for calling a fit good. Marsh, Hau and Wen (2004) argue against reading them as fixed rules, and with ordinal estimators the same indices tell a different story (Xia & Yang, 2019). They describe how far off a model is; the residuals show where.

Try it: what misfit does to fit indices

Six items, two factors correlated 0.3. The population also has things the model leaves out: item 3 loads on the second factor, and items 5 and 6 share an error. The widget fits the two-factor model by maximum likelihood.

Code

viewof fs_cross = Inputs.range([0, 0.5], {value: 0.3, step: 0.05, label: "Item 3's loading on factor 2 (not in the model)"})
viewof fs_cerr = Inputs.range([0, 0.3], {value: 0, step: 0.05, label: "Error covariance, items 5 and 6 (not in the model)"})
viewof fs_n = Inputs.range([100, 5000], {value: 500, step: 100, label: "Sample size N"})

Code

md`Expected χ² **${fs.idx.chisq.toFixed(1)}** on ${fs.idx.df} df (the 5% critical value is 15.5) · CFI **${fs.idx.cfi.toFixed(3)}** · RMSEA **${fs.idx.rmsea.toFixed(3)}** · SRMR **${fs.idx.srmr.toFixed(3)}** · estimated factor correlation **${fs.fit.Phi[0][1].toFixed(2)}** (true 0.3)`

At the start, the omitted cross-loading gives an RMSEA of 0.080 and a CFI of 0.944, and the factor correlation comes out at 0.48: the model routes item 3’s link to factor 2 through the correlation. The residuals point at item 3. Move N from 500 to 2,000: χ² rises from 34 to 111, and the other indices stay put. Now set the cross-loading to 0 and the error covariance to 0.2. CFI is 0.994 and RMSEA 0.026, inside the cutoffs, and the residual for items 5 and 6 is almost zero: the model raised their loadings to absorb the shared error, and the strain shows up between item 4 and items 1–3. A residual marks where the model strains, which isn’t always where the cause is.

Nested models can be compared directly: the one-factor model is the two-factor model with the factor correlation fixed at 1, so the difference in χ² tests that constraint. Robust estimators need a scaled difference (Satorra & Bentler, 2001), which lavTestLRT() computes. Prediction offers another view: Zhang et al. (2026) carry the IMV, which Model fit and out-of-sample prediction teaches, over to CFA with binary items.

Code

quiz({
  q: "A CFA fitted to 20,000 respondents has χ² = 450 on 100 degrees of freedom ($p < .001$) and CFI = 0.98. What should you make of the χ²?",
  options: [
    {text: "The model should be rejected", correct: false, why: "It is rejected as *exactly* true, but no model is exactly true, and with 20,000 respondents χ² detects tiny misfit. The misfit per respondent is small: RMSEA $= \\sqrt{(450 - 100)/(100 \\times 19{,}999)} \\approx 0.013$."},
    {text: "It says the model isn't exactly right, which with this N was always going to happen; how wrong, and where, is the question", correct: true, why: "χ² grows with $N$ for a fixed misfit (the widget shows it). RMSEA here is about 0.013 and CFI 0.98; the residuals show where the misfit is."},
    {text: "χ² is invalid with a sample this large", correct: false, why: "It is valid; it answers a question (exact fit) that large samples always answer 'no'."}
  ]
})

Omega: reliability from the model

Recall from Classical test theory and reliability: alpha is a lower bound on reliability, with equality under essential tau-equivalence. Factor analysis I read tau-equivalence as a one-factor model with equal loadings.

A factor model makes the bound unnecessary. Under a one-factor model with loadings \(\lambda_i\) and unique variances \(\psi_i\), the sum score has variance \((\sum_i \lambda_i)^2 + \sum_i \psi_i\), of which \((\sum_i\lambda_i)^2\) comes from the factor. Their ratio is omega (McDonald, 1999): \[\omega = \frac{(\sum_i \lambda_i)^2}{(\sum_i \lambda_i)^2 + \sum_i \psi_i}.\] When the loadings are equal, omega equals alpha (problem 1). When they aren’t, alpha is smaller, usually by a little (Revelle & Zinbarg, 2009).

Omega raises a sharper question when a scale has subscales. In a bifactor model (Holzinger & Swineford, 1937; Reise, 2012), every item loads on a general factor and on one specific factor for its subscale, with all factors uncorrelated. Omega total counts the general and specific factors as signal. Omega hierarchical, \(\omega_h\), counts only the general factor: it is the share of the total score’s variance that one common construct accounts for (Rodriguez, Reise & Haviland, 2016).

Try it: alpha, omega and omega hierarchical

Nine items in three subscales of three.

Code

viewof om_g = Inputs.range([0.3, 0.6], {value: 0.6, step: 0.05, label: "Average general-factor loading"})
viewof om_sp = Inputs.range([0, 0.3], {value: 0, step: 0.05, label: "Spread of general loadings (± within each subscale)"})
viewof om_s = Inputs.range([0, 0.4], {value: 0, step: 0.05, label: "Specific-factor loading"})

With equal loadings and no specific factors, all three are 0.835. Spread the loadings to ±0.2 and alpha falls a little below omega (0.833 against 0.841): the lower bound at work. Set the spread to 0 and the specific loading to 0.4: omega total rises to 0.886 and omega hierarchical falls to 0.771, because the subscales now carry reliable variance of their own. Alpha sits between the two.

Code

quiz({
  q: "A scale's bifactor model gives $\\omega_{\\text{total}} = 0.90$ and $\\omega_h = 0.80$. How is the variance of its total score divided?",
  options: [
    {text: "80% general factor, 10% specific factors, 10% unique", correct: true, why: "$\\omega_h$ is the general factor's share; $\\omega_{\\text{total}} - \\omega_h = 0.10$ is the specific factors' share; the remaining $1 - 0.90$ is unique variance."},
    {text: "90% general factor, 10% unique", correct: false, why: "That would make $\\omega_h = \\omega_{\\text{total}}$. The 0.10 gap between them belongs to the specific factors."},
    {text: "80% general factor, 20% unique", correct: false, why: "The specific factors are reliable variance too; omega total counts them, which is why it is 0.90."}
  ]
})

Item response theory as factor analysis

Recall from From the 1PL to the 4PL: the 2PL gives each item a slope \(a\) and a difficulty \(b\); mirt writes \(a\theta + d\), so \(b = -d/a\).

Factor analysis I treated an ordered response as a coarsened continuous one and factored polychoric correlations. Write that down for one item: a latent response \(y^* = \lambda\eta + \varepsilon\), with \(\text{var}(y^*) = 1\), and \(x = 1\) when \(y^*\) passes a threshold \(\tau\). Then \[\Pr(x = 1 \mid \eta) = \Phi\!\left(\frac{\lambda\eta - \tau}{\sqrt{1 - \lambda^2}}\right) = \Phi\big(a(\eta - b)\big), \quad a = \frac{\lambda}{\sqrt{1-\lambda^2}},\; b = \frac{\tau}{\lambda}.\] That is the normal-ogive 2PL. The factor model for dichotomous items is an item response model (Takane & de Leeuw, 1987), and with more categories each threshold becomes a category boundary of the graded response model (Samejima, 1969; Kamata & Bauer, 2008). For the logistic metric of mirt (Chalmers, 2012) and this course, multiply the slope by about 1.7 (Camilli, 1994, on where the constant comes from); this is the only place in the course it appears.

Try it: loading to slope

Code

md`Normal-ogive slope **${ls.irt.aNormal.toFixed(2)}**, logistic slope **${ls.irt.a.toFixed(2)}**, difficulty b = τ/λ = **${ls.irt.b.toFixed(2)}** (dashed line). Largest gap between the two curves: **${ls.maxDiff.toFixed(3)}**.`

At a loading of 0.7 and a threshold of 0.35, the slope is 0.98 (normal ogive) or 1.67 (logistic), the difficulty is 0.5, and the curves never differ by more than 0.01, wherever the sliders sit. Push the loading to 0.9 and the slope passes 3.5; at 0.95 it passes 5. Loadings stop at 1 and slopes don’t.

Go deeper: ordinal factor analysis is the graded response model

Assumptions: an item with ordered categories \(0, \dots, K\); a latent response \(y^* = \lambda^\top \eta + \varepsilon\) with \(\eta\) multivariate normal (variances 1), \(\varepsilon \sim N(0, 1 - h^2)\) independent of \(\eta\) and of the other items’ errors, where \(h^2 = \text{var}(\lambda^\top\eta)\) is the communality, so \(\text{var}(y^*) = 1\). The response is \(x = k\) when \(\tau_k < y^* \le \tau_{k+1}\), with \(\tau_0 = -\infty\) and \(\tau_{K+1} = \infty\).

Boundary curves. For \(k = 1, \dots, K\), \[\Pr(x \ge k \mid \eta) = \Pr(\varepsilon > \tau_k - \lambda^\top\eta) = \Phi\!\left(\frac{\lambda^\top\eta - \tau_k}{\sqrt{1 - h^2}}\right).\] With one factor this is \(\Phi\big(a(\eta - b_k)\big)\) with \(a = \lambda/\sqrt{1-\lambda^2}\) and \(b_k = \tau_k/\lambda\): one slope per item and one boundary per category, which is Samejima’s (1969) normal-ogive graded response model. Category probabilities are differences of adjacent boundary curves. With several factors, the slope on factor \(m\) is \(a_m = \lambda_m/\sqrt{1 - h^2}\).

Local independence. The errors are independent across items, so given \(\eta\) the responses are independent. Item response theory’s local independence is the factor model’s uncorrelated unique parts.

Logistic metric. \(\sup_z |\Phi(z) - \text{logit}^{-1}(1.702z)| < 0.01\) (the widget’s gap), so a logistic slope is about \(1.702a\); in mirt’s form, \(d_k = -1.702\,a\,b_k\).

Why the numbers differ in practice. The model is the same; the estimators are not. lavaan’s WLSMV estimates thresholds and polychoric correlations, which use the data one or two items at a time (Muthén, 1984; Flora & Curran, 2004). mirt maximizes the likelihood of whole response patterns. When the model holds they converge on the same values as the sample grows; in a finite sample they differ by sampling error, by the logistic approximation, and wherever the model misfits, since the two estimators weigh the misfit differently.

Choosing an estimator

Which estimator should a fit use? That depends on the categories and how they are used, the distance from normality, the sample size and what is missing. Three of lavaan’s estimators (tutorial) cover most scales.

ML (maximum likelihood), the default for numeric data, treats the items as continuous and multivariate normal (Li, 2016).
MLR keeps the ML estimates and makes the standard errors and the χ² robust to non-normality (Yuan & Bentler, 2000; Savalei, 2014, on how the corrections work).
WLSMV, lavaan’s choice when items are declared ordered, fits the ordinal model of the previous section to thresholds and polychoric correlations by diagonally weighted least squares (Muthén, 1984; lavaan tutorial, “Categorical data”). It assumes normal latent responses, not normal observed ones, and tolerates modest departures from that (Flora & Curran, 2004).

lavaan drops incomplete respondents by default; with ML or MLR, missing = "ml" gives full-information ML (FIML), which uses every observed response. With responses missing at random (Rubin, 1976), FIML was unbiased and more efficient than listwise or pairwise deletion in the simulations of Enders and Bandalos (2001). WLSMV has no FIML in lavaan; for ordered items, the full-information route is the same model as a graded response model in mirt (Bock & Aitkin, 1981).

When do they part ways? In simulations by Rhemtulla, Brosseau-Liard and Savalei (2012), the categorical estimator recovered loadings and standard errors better, especially with fewer than five categories, and robust ML did as well on factor correlations and fit; ML was the more sensitive to asymmetric thresholds, the categorical estimator to small samples and non-normal latent responses. With two or three categories, WLSMV’s χ² held its nominal rejection rate better (Beauducel & Herzberg, 2006). Li (2016) found WLSMV’s loadings less biased almost everywhere, but MLR’s standard errors better in small samples with non-normal latent responses. Summarized, with the caveat that each study speaks only for the conditions it simulated:

The responses	Where the simulations leaned	In `lavaan`
Continuous, roughly normal	ML	the default
Continuous but skewed, or five or more categories used fairly evenly	MLR	`estimator = "MLR"`
Fewer than five categories, or piled at one end	WLSMV	`ordered = TRUE`
Some responses missing	FIML with ML or MLR; with WLSMV, pairwise deletion or the graded model in `mirt`	`missing = "ml"` or `"pairwise"`
A few hundred respondents or fewer	MLR with five or more categories; else WLSMV, checked against MLR

The target matters too: omega is about the sum of the scored responses, so it needs a model of those responses even when the structure is fitted to polychorics. A second defensible estimator is a check on the first: where the two disagree, the choice matters for these data.

Code

quiz({
  q: "A three-category scale is answered by 800 respondents. Most choose the lowest category, and 5% of the responses are missing completely at random. Which statement about the options is right?",
  options: [
    {text: "ML on the complete cases makes no assumption these data could strain", correct: false, why: "ML treats three piled-up categories as normal numbers, which is where it was most sensitive in the simulations, and listwise deletion throws away every respondent with a single gap."},
    {text: "MLR with full-information ML uses every response, but still treats the three categories as numbers", correct: true, why: "FIML handles the missing responses, and MLR's standard errors are robust to non-normality, but the loadings come from treating the categories as numbers. With three categories piled at one end, the categorical estimator recovered loadings better in the simulations; fitting both shows whether that matters here."},
    {text: "WLSMV (items declared ordered) makes no distributional assumption", correct: false, why: "It assumes normal latent responses behind the categories, and in lavaan it has no FIML: with data missing completely at random, pairwise deletion is a reasonable companion, and the graded model in mirt uses every response."}
  ]
})

Simulate

The code simulates 1,000 respondents answering eight four-category items from two correlated factors, fits the CFA with lavaan and the graded model with mirt, and puts both on the logistic IRT scale. The first run takes 20–30 seconds.

With the seed as written, every slope from either route is within 0.3 of the truth, and the factor correlation comes out at 0.36 (true 0.4). Set np <- 300 and the slopes scatter; with 5,000 respondents they settle.

Download this code.

With real data

Grit: one factor or two?

Our main example is grit: the 12 items of the Grit-O (Duckworth, Peterson, Matthews & Kelly, 2007), answered online through Open Psychometrics. Its job is to test a designed two-facet structure. Download the code to run every chunk.

Code

library(lavaan)
# Each IRW table has a landing page (itemresponsewarehouse.org/tables/<name>/)
# with a CSV download. This link is pinned to one version of the data.
grit_url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v59_0.grit/rows?format=csv"
df <- read.csv(grit_url)
wide <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
grit_all <- as.data.frame(wide[, paste0("item_", 1:12)])
names(grit_all) <- paste0("g", 1:12)
grit_all <- grit_all[rowSums(!is.na(grit_all)) > 0, ]   # drop respondents with no responses
# Keying. The source (Open Psychometrics) coded 1 = "Very much like me" to
# 5 = "Not like me at all", and the IRW processing script reversed the six
# consistency-of-interest items (2, 3, 5, 7, 8, 11), which are worded against grit.
# So every item in the table runs toward LESS grit. Reversing all twelve makes
# higher = more grit, the same as reversing the perseverance items in the raw data.
grit_all <- 6 - grit_all
grit <- grit_all[complete.cases(grit_all), ]
per <- paste0("g", c(1, 4, 6, 9, 10, 12))   # perseverance of effort
con <- paste0("g", c(2, 3, 5, 7, 8, 11))    # consistency of interest
round(range(cor(grit)[lower.tri(diag(12))]), 2)   # no negative correlations left

[1] -0.03  0.56

On keying: the source coded 1 = “Very much like me” to 5 = “Not like me at all” (its codebook, in the data download), and the IRW’s processing script reversed the six consistency items, which are worded against grit. So every item in the table runs toward less grit, and reversing all twelve gives higher = more grit. No inter-item correlation is left below −0.03. Because the consistency items are exactly the ones worded against grit, facet and wording direction coincide here (problem 3).

Which estimator?

Code

# What the estimator choice depends on: the categories and how they are used,
# the shape of each item's distribution, and what is missing.
round(prop.table(table(unlist(grit))), 3)          # share of responses in each category


    1     2     3     4     5 
0.118 0.180 0.264 0.262 0.176

Code

skew <- function(x) mean((x - mean(x))^3) / sd(x)^3
round(range(sapply(grit, skew)), 2)                # item skewness, lowest and highest

[1] -0.91  0.26

Code

table(items_skipped = rowSums(is.na(grit_all)))    # respondents by number of items skipped

items_skipped
   0    1    2    3 
3103   48    4    3

Five categories, each holding 11.8% to 26.4% of the responses, and no item far from symmetric (skewness −0.91 to 0.26): where the simulations found the two estimators in agreement. So MLR, with WLSMV as the check. Only 55 of 3,158 respondents skipped any item, none more than three: we fit the 3,103 complete cases and check with FIML.

Recall from Factor analysis I: a one-factor model is what “the items are unidimensional” means, and it is a claim the data can reject.

Code

# lavaan syntax: "=~" reads "is measured by". Items not listed under a factor have
# a loading fixed at zero. By default lavaan fixes each factor's first loading to 1.
# MLR: maximum likelihood estimates, with standard errors and a test statistic
# that are robust to non-normality.
one <- "grit =~ g1 + g2 + g3 + g4 + g5 + g6 + g7 + g8 + g9 + g10 + g11 + g12"
two <- "per =~ g1 + g4 + g6 + g9 + g10 + g12
        con =~ g2 + g3 + g5 + g7 + g8 + g11"
f1 <- cfa(one, data = grit, estimator = "MLR")
f2 <- cfa(two, data = grit, estimator = "MLR")
idx <- c("chisq.scaled", "df", "cfi.robust", "tli.robust", "rmsea.robust", "srmr")
round(rbind(one_factor = fitMeasures(f1, idx), two_factors = fitMeasures(f2, idx)), 3)

            chisq.scaled df cfi.robust tli.robust rmsea.robust  srmr
one_factor      2418.471 54      0.758      0.704        0.127 0.095
two_factors     1052.292 53      0.898      0.873        0.084 0.066

Code

lavTestLRT(f1, f2)   # scaled difference: the one-factor model is the two-factor model with the correlation fixed at 1

	Df	AIC	BIC	Chisq	Chisq diff	RMSEA	Df diff	Pr(>Chisq)
f2	53	106532.0	106683	1211.399	NA	NA	NA	NA
f1	54	108098.1	108243	2779.444	1498.801	0.7106299	1	0

Code

lavInspect(f2, "cor.lv")

      per   con
per 1.000      
con 0.612 1.000

The data reject it: one factor gives a robust CFI of 0.758 and an RMSEA of 0.127. The designed facets improve both (0.898 and 0.084, SRMR 0.066), and the scaled difference test is decisive, though by the conventional cutoffs two factors still fall short. The facets correlate 0.61.

Code

# The same two-factor model three ways: plain ML and MLR on the complete cases
# (same estimates, different test statistic), and MLR with full-information ML
# (missing = "ml") on every respondent who answered at least one item.
f2_ml <- cfa(two, data = grit)
f2_fiml <- cfa(two, data = grit_all, estimator = "MLR", missing = "ml")
both <- function(f, i) setNames(fitMeasures(f, i), c("cfi", "rmsea"))
round(rbind(ML = c(n = nobs(f2_ml), both(f2_ml, c("cfi", "rmsea")), factor_cor = lavInspect(f2_ml, "cor.lv")[1, 2]),
            MLR = c(n = nobs(f2), both(f2, c("cfi.robust", "rmsea.robust")), factor_cor = lavInspect(f2, "cor.lv")[1, 2]),
            MLR_FIML = c(n = nobs(f2_fiml), both(f2_fiml, c("cfi.robust", "rmsea.robust")),
                         factor_cor = lavInspect(f2_fiml, "cor.lv")[1, 2])), 3)

            n   cfi rmsea factor_cor
ML       3103 0.897 0.084      0.612
MLR      3103 0.898 0.084      0.612
MLR_FIML 3158 0.898 0.084      0.616

Code

round(max(abs(coef(f2)[1:12] - coef(f2_fiml)[1:12])), 3)   # largest change in a loading with FIML

[1] 0.01

Plain ML gives the same estimates and nearly the same indices; FIML on all 3,158 respondents moves no loading by more than 0.01.

Code

# Two ways to set the scale of each factor: fix the first loading to 1 (lavaan's
# default) or fix the factor's variance to 1 (std.lv = TRUE). Same fit, different
# numbers for the loadings.
f2_var <- cfa(two, data = grit, estimator = "MLR", std.lv = TRUE)
rbind(first_loading_1 = fitMeasures(f2, c("chisq", "df")),
      variance_1      = fitMeasures(f2_var, c("chisq", "df")))

                   chisq df
first_loading_1 1211.399 53
variance_1      1211.399 53

Code

round(cbind(first_loading_1 = coef(f2)[1:6], variance_1 = coef(f2_var)[1:6]), 2)

         first_loading_1 variance_1
per=~g4             0.80       0.55
per=~g6             1.33       0.44
per=~g9             1.40       0.73
per=~g10            1.56       0.77
per=~g12            1.28       0.86
con=~g3             1.14       0.70

Both ways of setting the scale give the same χ² on 53 df; only the loadings’ units change.

Code

# Omega for each facet from the two-factor model: (sum of loadings)^2 over the
# model-implied variance of the facet's sum score. Alpha from its usual formula.
est <- lavInspect(f2_var, "est")
omega <- function(items, factor) {
  l <- est$lambda[items, factor]; th <- est$theta[items, items]
  sum(l)^2 / (sum(l)^2 + sum(th))
}
alpha <- function(x) { k <- ncol(x); k / (k - 1) * (1 - sum(apply(x, 2, var)) / var(rowSums(x))) }
round(rbind(omega = c(perseverance = omega(per, "per"), consistency = omega(con, "con")),
            alpha = c(alpha(grit[per]), alpha(grit[con]))), 3)

      perseverance consistency
omega        0.777       0.811
alpha        0.772       0.809

Code

round(est$lambda[per, "per"], 2)   # unequal loadings: the items aren't tau-equivalent

  g1   g4   g6   g9  g10  g12 
0.55 0.44 0.73 0.77 0.86 0.70

Omega and alpha nearly agree (0.777 and 0.772 for perseverance): the loadings run from 0.44 to 0.86, but the gap is small, as in the widget.

Now the items as ordered categories, and the same model as a graded response model in mirt:

Code

# The same two-factor model with the items declared ordered: lavaan then fits the
# polychoric correlations (estimator WLSMV). Loadings are for the latent responses.
f2_ord <- cfa(two, data = grit, ordered = TRUE, std.lv = TRUE)
idx_s <- c("chisq.scaled", "df.scaled", "cfi.scaled", "rmsea.scaled", "srmr")
round(fitMeasures(f2_ord, idx_s), 3)

chisq.scaled    df.scaled   cfi.scaled rmsea.scaled         srmr 
    2131.609       53.000        0.906        0.112        0.070

Code

lam <- lavInspect(f2_ord, "std")$lambda
lam <- rowSums(lam)[c(per, con)]          # each item's one nonzero loading
tau <- matrix(lavInspect(f2_ord, "est")$tau, ncol = 4, byrow = TRUE,
              dimnames = list(rownames(lavInspect(f2_ord, "est")$lambda), NULL))[c(per, con), ]
# Convert: normal-ogive slope lambda / sqrt(1 - lambda^2), times 1.702 for the
# logistic metric; difficulty (category boundary) b = tau / lambda.
a_from_cfa <- 1.702 * lam / sqrt(1 - lam^2)
b_from_cfa <- tau / lam

Code

library(mirt)
spec <- mirt.model("per = 1, 4, 6, 9, 10, 12
                    con = 2, 3, 5, 7, 8, 11
                    COV = per*con")
grm <- mirt(grit, spec, itemtype = "graded", verbose = FALSE)
cf <- coef(grm, simplify = TRUE)
a_mirt <- rowSums(cf$items[, c("a1", "a2")])[c(per, con)]
# mirt writes each boundary as a*theta + d; the difficulty is b = -d/a.
b_mirt <- -cf$items[c(per, con), c("d1", "d2", "d3", "d4")] / a_mirt
round(cbind(loading = lam, a_from_cfa, a_mirt), 2)

    loading a_from_cfa a_mirt
g1     0.55       1.13   1.32
g4     0.40       0.74   0.81
g6     0.74       1.87   2.05
g9     0.78       2.12   1.79
g10    0.69       1.63   1.61
g12    0.71       1.73   1.90
g2     0.60       1.28   1.36
g3     0.61       1.33   1.41
g5     0.73       1.83   1.95
g7     0.78       2.13   2.37
g8     0.81       2.36   2.09
g11    0.46       0.87   1.00

Code

round(c(cor_a = cor(a_from_cfa, a_mirt), cor_b = cor(as.vector(b_from_cfa), as.vector(b_mirt)),
        largest_b_gap = max(abs(b_from_cfa - b_mirt)),
        factor_cor_lavaan = lavInspect(f2_ord, "cor.lv")[1, 2],
        factor_cor_mirt = cf$cov[1, 2]), 2)

            cor_a             cor_b     largest_b_gap factor_cor_lavaan 
             0.94              1.00              0.29              0.62 
  factor_cor_mirt 
             0.61

WLSMV tells the same story in different numbers (CFI 0.906, RMSEA 0.112). The converted slopes track mirt’s (correlation 0.94; single items differ by up to 0.33), the difficulties agree closely (correlation 1.00), and the factor correlation is 0.62 by one route and 0.61 by the other. Item 11 has the weakest loading, 0.46; Dimensionality and multidimensional IRT meets it again.

The DASS-21: three scales or one?

Our contrast table is neurodegenerative_huizinga_2019_dass: the Dutch DASS-21, answered by a community sample of adults recruited by Huizinga et al. (2020) for a study of visual complaints. Its job: a structure whose designed scales may not be distinct. The DASS-21 scores depression, anxiety and stress with seven items each (Lovibond & Lovibond, 1995; Henry & Crawford, 2005). Every item is worded toward more distress, and the scales follow the DASS scoring template.

Code

dass_url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v59_0.neurodegenerative_huizinga_2019_dass/rows?format=csv"
df <- read.csv(dass_url)
df$resp <- suppressWarnings(as.numeric(df$resp))   # missing responses arrive as text
wide <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
dass <- as.data.frame(wide[, paste0("DASS.", 1:21)])
names(dass) <- paste0("d", 1:21)
table(items_skipped = rowSums(is.na(dass)))   # all or nothing: no partial missingness

items_skipped
   0   21 
1246  215

Code

dass <- dass[complete.cases(dass), ]
# 0 = did not apply to me at all ... 3 = applied to me very much. Every item is
# worded toward more distress, so no keying is needed.
table(response = unlist(dass))

response
    0     1     2     3 
20031  5103   851   181

Code

round(range(sapply(dass, skew)), 2)            # item skewness, lowest and highest

[1] 1.24 3.38

Code

# Scales from the DASS scoring template (Lovibond & Lovibond, 1995).
dep <- paste0("d", c(3, 5, 10, 13, 16, 17, 21))
anx <- paste0("d", c(2, 4, 7, 9, 15, 19, 20))
stress <- paste0("d", c(1, 6, 8, 11, 12, 14, 18))
m3 <- paste("dep =~", paste(dep, collapse = " + "), "\n",
            "anx =~", paste(anx, collapse = " + "), "\n",
            "str =~", paste(stress, collapse = " + "))
m1 <- paste("distress =~", paste(c(dep, anx, stress), collapse = " + "))
mbi <- paste(m1, "\n", m3)   # bifactor: every item on the general factor and on its scale

Which estimator this time? The 215 respondents without DASS responses skipped all 21 items and the rest answered all 21, so FIML has nothing to recover: we keep 1,246 complete cases. Four categories, most responses 0, only 181 of 26,166 at 3, every item skewed right (skewness 1.24 to 3.38): the conditions where the simulations leaned toward WLSMV. We use it, except for omega (MLR).

Code

predict({
  q: "The DASS-21 is scored as three scales, and a three-factor model fits these data well (CFI 0.975). How strongly will the three factors correlate?",
  choices: ["About 0.3", "About 0.5", "About 0.7", "Above 0.8"],
  answer: "Above 0.8",
  reveal: "Between 0.84 and 0.87. In this sample the three scales are close to one construct: a respondent high on one tends to be high on all three."
})

Code

# Four categories with most responses at 0: fit the polychoric correlations.
d3 <- cfa(m3, data = dass, ordered = TRUE, std.lv = TRUE)
round(lavInspect(d3, "cor.lv"), 2)

     dep  anx  str
dep 1.00          
anx 0.84 1.00     
str 0.86 0.87 1.00

Code

d1 <- cfa(m1, data = dass, ordered = TRUE, std.lv = TRUE)
dbi <- cfa(mbi, data = dass, ordered = TRUE, std.lv = TRUE, orthogonal = TRUE)
round(rbind(one_factor = fitMeasures(d1, idx_s), three_factors = fitMeasures(d3, idx_s),
            bifactor = fitMeasures(dbi, idx_s)), 3)

              chisq.scaled df.scaled cfi.scaled rmsea.scaled  srmr
one_factor        1376.362       189      0.948        0.071 0.060
three_factors      744.101       186      0.975        0.049 0.043
bifactor           437.745       168      0.988        0.036 0.032

The factors correlate 0.84 to 0.87. One factor fits a little worse (CFI 0.948 against 0.975), the bifactor model best (0.988). How much of the total score is one general distress factor?

Code

# Omega for the sum score needs a model of the scored responses, so here the items
# are treated as numbers (MLR: maximum likelihood with robust standard errors).
bi <- cfa(mbi, data = dass, estimator = "MLR", std.lv = TRUE, orthogonal = TRUE)
L <- lavInspect(bi, "est")$lambda
th <- diag(lavInspect(bi, "est")$theta)
total_var <- sum(colSums(L)^2) + sum(th)
omega_h_sub <- c(dep = sum(L[dep, "dep"])^2 / (sum(colSums(L[dep, ])^2) + sum(th[dep])),
                 anx = sum(L[anx, "anx"])^2 / (sum(colSums(L[anx, ])^2) + sum(th[anx])),
                 str = sum(L[stress, "str"])^2 / (sum(colSums(L[stress, ])^2) + sum(th[stress])))
round(c(omega_total = sum(colSums(L)^2) / total_var,
        omega_hierarchical = sum(L[, "distress"])^2 / total_var,
        alpha = alpha(dass)), 3)

       omega_total omega_hierarchical              alpha 
             0.944              0.879              0.933

Code

round(omega_h_sub, 3)   # share of each scale's sum-score variance that is specific to it

  dep   anx   str 
0.240 0.230 0.002

Omega total is 0.944 and omega hierarchical 0.879: one general factor accounts for almost all of the total score’s reliable variance. Within each scale, the specific factor’s share of the scale score is 0.24 for depression, 0.23 for anxiety and 0.002 for stress; the stress items carry almost nothing beyond general distress.

Should a clinic report three DASS scores or one? The subscales add information only through their specific factors, which are small here, and the DASS authors describe the total as a sensible composite (DASS FAQ, no. 26). In the case of a community sample like this one, the total may be sufficient for some purposes, with the depression and anxiety scores reported only alongside it. A clinical sample could give a different answer. Would women and men with the same level of distress answer the DASS alike? Measurement invariance takes this table up there.

Problems

Omega and alpha. Show that under a one-factor model with equal loadings \(\lambda\) and equal unique variances, omega equals alpha. Then show that with unequal loadings (and uncorrelated errors) alpha is smaller.
The probit link. Starting from \(y^* = \lambda\eta + \varepsilon\) with \(\text{var}(y^*) = 1\), derive the normal-ogive slope and difficulty of a dichotomized item. Use your result on grit: which facet’s items discriminate better, and does mirt agree?
A patch or a finding? (a) Run modindices() on the two-factor grit model; among the largest suggestions is a correlated error between items 3 and 11. Add it. How much does fit improve? Given the codebook’s stems (“My interests…”, “I become…”), is it a finding or a patch? Can any model of these data separate consistency of interest from wording direction (From construct map to items)? (b) Dichotomize the keyed items (4 or 5 against the rest), fit the one- and two-factor binary CFAs, and compare them by the IMV (Model fit and out-of-sample prediction) on held-out respondents with five-fold cross-validation, as Zhang et al. (2026) do. Does prediction agree with the fit indices, item by item?
Two estimators, one DASS. Refit the three-factor DASS model with MLR and compare its standardized loadings and factor correlations with the WLSMV fit. Which move, and which way? Is that what Rhemtulla et al. (2012) and Li (2016) would lead you to expect?
Design. A new three-category scale (never, sometimes, often) will go to about 250 respondents; a pilot suggests most answers will be “never” and 1 in 10 respondents will skip an item. Which estimator and treatment of missing responses would you plan on, what would you fit as a check, and what would change your mind?
Challenge (open). Credé, Tynan and Harms (2017) argue from a meta-analysis that grit overlaps heavily with conscientiousness. What confirmatory model would test that, and what data would you need? I don’t know of a settled answer.

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

Jöreskog, K. G. (1969). A general approach to confirmatory maximum likelihood factor analysis. Psychometrika, 34(2), 183–202. https://doi.org/10.1007/bf02289343. Where CFA starts.
Rosseel, Y. (2012). lavaan: An R package for structural equation modeling. Journal of Statistical Software, 48(2). https://doi.org/10.18637/jss.v048.i02. The package and its syntax; Structural equation modeling with lavaan adds paths between factors.
Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling, 6(1), 1–55. https://doi.org/10.1080/10705519909540118. Where the cutoffs come from, and under which conditions.
Marsh, H. W., Hau, K.-T., & Wen, Z. (2004). In search of golden rules. Structural Equation Modeling, 11(3), 320–341. https://doi.org/10.1207/s15328007sem1103_2. Argues against reading the cutoffs as rules.
Zhang, L., Rahal, C., Kanopka, K., Ulitzsch, E., Zhang, Z., & Domingue, B. W. (2026). Evaluating model predictive performance in confirmatory factor analysis with binary outcomes using the InterModel Vigorish. Multivariate Behavioral Research. https://doi.org/10.1080/00273171.2026.2645212. Prediction on held-out data as a complement to fit indices for CFA with binary items, down to the single item (problem 3).
Xia, Y., & Yang, Y. (2019). RMSEA, CFI, and TLI in structural equation modeling with ordered categorical data: The story they tell depends on the estimation methods. Behavior Research Methods, 51(1), 409–428. https://doi.org/10.3758/s13428-018-1055-2. Why the grit indices changed with the estimator.
McDonald, R. P. (1999). Test theory: A unified treatment. Erlbaum. https://doi.org/10.4324/9781410601087. Omega, and test theory built on the factor model.
Revelle, W., & Zinbarg, R. E. (2009). Coefficients alpha, beta, omega, and the glb: Comments on Sijtsma. Psychometrika, 74(1), 145–154. https://doi.org/10.1007/s11336-008-9102-z. Omega hierarchical beside the other coefficients.
Reise, S. P. (2012). The rediscovery of bifactor measurement models. Multivariate Behavioral Research, 47(5), 667–696. https://doi.org/10.1080/00273171.2012.715555. What a bifactor model can say about subscales.
Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models: Calculating and interpreting statistical indices. Psychological Methods, 21(2), 137–150. https://doi.org/10.1037/met0000045. Omega hierarchical and its relatives, computed and read.
Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using structural equation modeling: An alternative to coefficient alpha. Psychometrika, 74(1), 155–167. https://doi.org/10.1007/s11336-008-9099-3. Omega for ordinal items, on the scale of the sum score.
Takane, Y., & de Leeuw, J. (1987). On the relationship between item response theory and factor analysis of discretized variables. Psychometrika, 52(3), 393–408. https://doi.org/10.1007/bf02294363. The equivalence, proved.
Kamata, A., & Bauer, D. J. (2008). A note on the relation between factor analytic and item response theory models. Structural Equation Modeling, 15(1), 136–153. https://doi.org/10.1080/10705510701758406. The conversions under each parameterization.
Rhemtulla, M., Brosseau-Liard, P. É., & Savalei, V. (2012). When can categorical variables be treated as continuous? A comparison of robust continuous and categorical SEM estimation methods under suboptimal conditions. Psychological Methods, 17(3), 354–373. https://doi.org/10.1037/a0029315. The five-category guideline.
Flora, D. B., & Curran, P. J. (2004). An empirical evaluation of alternative methods of estimation for confirmatory factor analysis with ordinal data. Psychological Methods, 9(4), 466–491. https://doi.org/10.1037/1082-989X.9.4.466. How WLSMV behaves.
Li, C.-H. (2016). Confirmatory factor analysis with ordinal data: Comparing robust maximum likelihood and diagonally weighted least squares. Behavior Research Methods, 48(3), 936–949. https://doi.org/10.3758/s13428-015-0619-7. MLR against WLSMV across numbers of categories, latent distributions and sample sizes.
Enders, C. K., & Bandalos, D. L. (2001). The relative performance of full information maximum likelihood estimation for missing data in structural equation models. Structural Equation Modeling, 8(3), 430–457. https://doi.org/10.1207/S15328007SEM0803_5. FIML against listwise and pairwise deletion.
The lavaan tutorial pages on estimators and categorical data: what each estimator and missing option does.
Measurement invariance, the next lesson: the same factor models fitted to two groups at once.
From the IRW: the CFA walkthrough with lavaan, and the dimensionality vignette on how often one factor holds.

Data sources

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

grit: Duckworth, A. L., Peterson, C., Matthews, M. D., & Kelly, D. R. (2007). Grit: perseverance and passion for long-term goals. Journal of personality and social psychology, 92(6), 1087. https://doi.org/10.1037/0022-3514.92.6.1087. Licence: CC BY 4.0.
neurodegenerative_huizinga_2019_dass: Huizinga, F., Heutink, J., de Haan, G. A., van der Lijn, I., van der Feen, F. E., Vrijling, A. C., … & Koerts, J. (2020). The development of the Screening of Visual Complaints questionnaire for patients with neurodegenerative disorders: Evaluation of psychometric features in a community sample. PLoS One, 15(4), e0232232. https://doi.org/10.1371/journal.pone.0232232. Licence: CC0 1.0.

For instructors

Source: new material. EDUC 252 did not cover confirmatory factor analysis; the lesson was added after the landscape analysis (notes/landscape-2026-09-24.md), where CFA appeared in 6 of 11 first-course syllabi. Multigroup invariance moved here from Differential item functioning (Ben, 09-24), and on to its own lesson, Measurement invariance, when this one was split (Ben, 09-25, #58). The estimator section was added at the same time, at a colleague’s suggestion: the choice among ML, MLR and WLSMV is made explicitly for both tables.
Timing: about 60 minutes for the core ideas and widgets; the real-data section is a 30-minute lab (the two-dimensional graded model takes under a minute locally).
Placement: after From the 1PL to the 4PL, although it sits in the factor-analysis module, because the IRT connection needs the 2PL. lavaan syntax is introduced here; Structural equation modeling with lavaan builds on it.
Held back: multigroup CFA and measurement invariance (Measurement invariance); paths between latent variables and multigroup SEM (Structural equation modeling with lavaan); multidimensional IRT (Dimensionality and multidimensional IRT); the graded response model in its own right (Models for polytomous responses); invariance item by item with a matching variable (Differential item functioning); omega for reported scores (What does a score mean?).
Data: the DASS is in the public domain (DASS FAQ, no. 3); no item is quoted here. The Grit-O items are described, not reproduced. Grit is fitted with MLR (the FIML and plain-ML fits are shown as checks) and with WLSMV; the DASS structure with WLSMV, and its omega with MLR.
Solutions: held back in solutions/fa-confirmatory.qmd (not published).