Structural equation modeling with lavaan
What this is for
Does grit go with feeling good? A regression of one sum score on another answers a smaller question than it seems to: it relates two error-laden measures, and error shrinks and shuffles the coefficients. Factor analysis II modelled each construct as a factor behind its items. A structural equation model (SEM) keeps those factors and adds regressions between them, so the relations are estimated between the constructs rather than between their scores (Jöreskog, 1970; Bollen, 1989). That buys a correction for measurement error, and it costs assumptions: a measurement model that must be right before the paths mean anything, and paths that look causal whether or not they are. We use the Grit-O and the PANAS, answered by the same children and adolescents in the Florida Twin Project.
Goals
By the end of this lesson you should be able to:
- say what the measurement and structural parts of an SEM are, and write both in
lavaansyntax; - explain how measurement error biases regressions on sum scores, including with two correlated predictors, and why a latent regression doesn’t carry that bias when its model holds;
- check by counting whether a structural model can be identified, and say which part of the model a fit index describes;
- fit and inspect a measurement model before the structural one, and compare latent paths with sum-score regressions on real data;
- say why a well-fitting path model is not evidence for its causal direction, and write down an equivalent model.
Core ideas
A measurement model and a structural model
An SEM has two parts. The measurement model is a confirmatory factor model: each factor \(\eta\) is measured by its items, \(x_i = \nu_i + \lambda_i \eta + \varepsilon_i\). The structural model relates the factors to each other with regressions, \[\eta_{\text{outcome}} = \beta_1 \eta_1 + \beta_2 \eta_2 + \zeta,\] where \(\zeta\) is the part of the outcome the predictors leave unexplained. Path analysis, regressions among observed variables drawn as arrows, goes back to Wright (1934); Jöreskog (1970) joined it to the factor model and fitted both by maximum likelihood.
Our model has four factors: two facets of grit (perseverance of effort, consistency of interest) and two kinds of affect (positive, negative). The question is how the grit facets relate to affect, with each construct measured by its items.
Why a latent path differs from a regression on sum scores
With one predictor, the standardized regression slope is the correlation, so it shrinks by the same factor. Error in the outcome and error in the predictor act differently on the unstandardized slope, though: error in \(Y\) adds noise around the line without tilting it, while error in \(X\) flattens it by the factor \(\rho_{XX'}\) (problem 1).
With two correlated predictors the damage is no longer a simple shrinkage. Controlling for a predictor measured with error controls for it only in part, so what is left of it leaks into the other predictor’s coefficient (Cole & Preacher, 2014; Westfall & Yarkoni, 2016). An SEM avoids both problems when its model holds: the paths run between the factors, which carry no measurement error, and the correction for attenuation is built into the model rather than applied by hand. A regression on sum scores rests on different assumptions: that the scores are reliable enough for the shrinkage and the leak to be small, or that the score itself, not the construct, is what will be used, as with a cut score.
Try it: two predictors, measured with error
Set the true standardized paths from two latent predictors to a latent outcome, and the reliability of the sum score that measures each. The bars compare the true paths with the standardized coefficients of the regression on sum scores.
At the start, predictor 2 has no effect at all, yet its sum score gets a coefficient of 0.07, and predictor 1’s path of 0.5 comes out at 0.39. Lower predictor 1’s reliability to 0.5 and predictor 2’s spurious coefficient doubles, to 0.14: the less well predictor 1 is measured, the more of it predictor 2 absorbs. Set the correlation to 0 and the leak stops, leaving plain attenuation. Set every reliability to 1 and the bars agree whatever the paths.
Specifying the model in lavaan
Two more operators finish the language: ~ reads “is regressed on”, and ~~ declares a variance or covariance. The structural part of our model is two lines added under the four =~ lines:
model <- "per =~ g1 + g4 + g6 + g9 + g10 + g12
con =~ g2 + g3 + g5 + g7 + g8 + g11
pa =~ p1 + p3 + p5 + p9 + p10 + p12 + p14 + p16 + p17 + p19
na =~ p2 + p4 + p6 + p7 + p8 + p11 + p13 + p15 + p18 + p20
pa ~ per + con
na ~ per + con"
fit <- sem(model, data = d, estimator = "MLR")sem() lets the two predictors correlate and, by default, lets the two outcomes’ residuals \(\zeta\) correlate too.
Can the structural part be identified? The measurement model already estimates the six correlations among four factors. The structural model re-expresses them with the predictors’ correlation, the paths and any residual covariance. With more than six such parameters, some can’t be recovered from the data. Counting is necessary but not sufficient for identification (Bollen, 1989).
The estimator matters as it did in CFA. Maximum likelihood assumes multivariate normal items. MLR keeps the ML estimates and makes the standard errors and the test statistic robust to non-normality (Yuan & Bentler, 2000); WLSMV treats items as ordered categories and suits items with few or skewed categories (Rhemtulla, Brosseau-Liard & Savalei, 2012).
Try it: build the structural model
Choose the paths. The widget draws them, writes the lavaan lines (the measurement lines above stay the same) and counts the structural degrees of freedom.
Dashed curves are correlations: grey between the predictors, orange between the outcomes’ residuals. With all four paths and the residual covariance, the count is 0: the structural model is only a new way of writing the six factor correlations, and it fits exactly as well as the measurement model alone. Drop “pa ~ con” and a structural claim appears that the data can test. Add “na ~ pa” with everything else on and the count goes to −1: the direct path and the residual covariance can’t both be recovered.
Fit, and fit of which part
A global fit index compares the whole implied covariance matrix with the data. In an SEM most of that matrix is the measurement model: in ours, 32 items supply 528 variances and covariances, and the structural claims concern at most six factor correlations. Misfit among the items can swamp the structure, and a structural part with few degrees of freedom of its own can be wrong without the global indices showing it. Anderson and Gerbing (1988) proposed two steps: fit the measurement model alone, with every factor correlated, and settle it; then fit the structural model and test its constraints against the measurement model. The cutoffs of Hu and Bentler (1999) come from simulations of factor models, and a structural model inherits whatever misfit its measurement model carries.
Fitting everything in one step gives the same estimates. In the case that the measurement model is well established in similar samples, one step may be sufficient for some purposes; what it can’t do is say which part a misfit belongs to.
Fitting and inspecting the measurement model first says how well the items measure each factor before any structural path is read. A path between factors that the items don’t measure well is an estimate of something, but not of the relation the diagram names, and the structural fit indices won’t say so.
Factors or composites? Everything here models a construct as a common factor, the source of what its items share. Composite-based SEM (partial least squares, consistent PLS, generalized structured component analysis) models a construct as a weighted sum of its items instead. It is a different model with different assumptions, closer to a construct that is defined by its indicators, as an index of socioeconomic status is. In R it is fitted by cSEM (Rademaker & Schuberth, 2020; on the model, Schuberth, Henseler & Dijkstra, 2018), and in jamovi by the Composite-SEM module, a small project without a licence file as of September 2026. For point-and-click factor-based SEM, jamovi’s SEMLj module builds lavaan models from menus.
What a path can’t tell you
An arrow in a path diagram is a regression coefficient. Whether it is a causal effect depends on assumptions the fit doesn’t test: no omitted common causes, the right direction, and a design that supports both (Bollen & Pearl, 2013). And for most models there are equivalent models, with different arrows, that imply exactly the same covariance matrix and so fit identically on any data (MacCallum, Wegener, Uchino & Fabrigar, 1993). Whether the constructs themselves cause their indicators is a further question, which Validity as causation takes up.
Try it: reverse an arrow
Three observed, standardized variables with the correlations you set, and 500 respondents. Each model leaves out one link and so makes one testable claim (1 degree of freedom).
The chain, the reversed chain and the fork all give χ² = 8.0 on 1 df, and they stay equal wherever the sliders sit: each says only that X and Y are unrelated given M. Set the X–Y correlation to 0.20, the product of the other two, and all three fit perfectly. They tell three causal stories (M a mediator one way, M a mediator the other way, M a common cause) that no amount of data of this kind can separate. The collider claims instead that X and Y are unrelated, and its χ² of 47.2 differs.
Simulate
The code simulates a latent predictor and a latent outcome with a standardized path of 0.5, each measured by five items, at three levels of item quality. For each it compares the correlation of the sum scores (the standardized slope of a one-predictor regression), Spearman’s prediction of it, and the path from sem(). The first run takes 20–30 seconds while lavaan loads.
With the seed as written, the sum-score slope falls from 0.45 to 0.24 as the loadings drop from 0.8 to 0.4, within 0.01 of Spearman’s prediction each time, while the SEM path stays between 0.50 and 0.56. The SEM estimate is unbiased, not free: with loadings of 0.4 it moves more from sample to sample. Change the seed, or set np <- 200, to see how much; with np <- 5000 (locally) it settles near 0.5.
With real data
Two tables from one study
The Florida Twin Project on Reading, Behavior, and Environment surveyed twins in Florida over three waves (Hart, Schatschneider & Taylor, 2021; on the project, Hart et al., 2013). Our tables are florida_twins_grit, the 12-item Grit-O (Duckworth, Peterson, Matthews & Kelly, 2007), and florida_twins_panas, the 20-item PANAS (Watson, Clark & Tellegen, 1988): ten adjectives for positive affect and ten for negative, each rated from 1 (very slightly or not at all) to 5 (extremely) for how the respondent feels in general. They are one data source, joined by respondent, with one job: grit supplies the predictors and affect the outcomes of the structural model. 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. These links are pinned to one version of the data.
base <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v60_0."
grit_long <- read.csv(paste0(base, "florida_twins_grit/rows?format=csv"))
panas_long <- read.csv(paste0(base, "florida_twins_panas/rows?format=csv"))
table(grit = grit_long$wave[!duplicated(grit_long$id)]) # grit was asked in wave 2 onlygrit
2
782 Code
panas_long <- panas_long[panas_long$wave == 2, ] # so we use PANAS from wave 2
wide <- function(df, items) {
w <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
as.data.frame(w[, items])
}
grit <- wide(grit_long, paste0("qgrit", 1:12)); names(grit) <- paste0("g", 1:12)
panas <- wide(panas_long, paste0("panas", 1:20, "t")); names(panas) <- paste0("p", 1:20)
# Keying. Grit runs 1 = "Very much like me" to 5 = "Not like me at all" (the IRW's
# option text). The consistency-of-interest items are worded against grit ("My
# interests change..."), so "not like me" already means more grit; the six
# perseverance items are reversed. Higher = more grit on every item.
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
grit[per] <- 6 - grit[per]
# PANAS: 1 = "very slightly or not at all" to 5 = "extremely"; no keying needed.
# Scales as in Watson, Clark & Tellegen (1988).
pa <- paste0("p", c(1, 3, 5, 9, 10, 12, 14, 16, 17, 19)) # positive affect
na <- paste0("p", c(2, 4, 6, 7, 8, 11, 13, 15, 18, 20)) # negative affect
d <- merge(grit, panas, by = "row.names")
d$id <- as.numeric(d$Row.names)
# Twins: IRW ids come in pairs ending 00 and 01 (31700, 31701), so id %/% 100 is
# the family. Inferred from the id pattern, not from a codebook.
d$family <- d$id %/% 100
d$age <- tapply(grit_long$cov_age, grit_long$id, function(x) x[1])[as.character(d$id)]
c(both_tables = nrow(d), complete = sum(complete.cases(d[c(per, con, pa, na)])))both_tables complete
778 618 Code
d <- d[complete.cases(d[c(per, con, pa, na)]), ]
table(twins_in_family = table(d$family))twins_in_family
1 2
110 254 Code
r <- cor(d[c(per, con)])
round(r[per, c("g3", "g11")], 2) # the two items that don't line up with the rest g3 g11
g1 -0.06 -0.21
g4 -0.10 -0.07
g6 -0.05 -0.18
g9 -0.02 -0.19
g10 -0.01 -0.33
g12 0.00 -0.25Grit was asked in wave 2 only, so we use the wave-2 PANAS. The grit items run from 1 = “Very much like me” to 5 = “Not like me at all”. The six consistency items are worded against grit, so their high end already means more grit, and we reverse the six perseverance items. (The Open Psychometrics grit table in Factor analysis II arrives with its consistency items already reversed, so there all twelve are reversed; the result, higher = more grit, is the same.) After keying, two consistency items sit apart: item 11, on taking up new interests every few months, correlates negatively with every perseverance item, and item 3 barely correlates with them (−0.10 to 0.00). Dimensionality and multidimensional IRT meets item 11 again.
Twins share a family, and the IRW ids come in pairs (31700 and 31701), so we take id %/% 100 as the family. That is read off the pattern of the ids, not a codebook.
Step 1: the measurement model
Code
# Step 1: the measurement model alone. Four factors, each item on one, all four
# factors correlated. MLR: maximum likelihood with standard errors and a test
# statistic that are robust to non-normal responses; cluster = "family" makes them
# robust to twins' responses being correlated as well.
mm <- paste0("per =~ ", paste(per, collapse = " + "), "\n",
"con =~ ", paste(con, collapse = " + "), "\n",
"pa =~ ", paste(pa, collapse = " + "), "\n",
"na =~ ", paste(na, collapse = " + "))
cat(mm)per =~ g1 + g4 + g6 + g9 + g10 + g12
con =~ g2 + g3 + g5 + g7 + g8 + g11
pa =~ p1 + p3 + p5 + p9 + p10 + p12 + p14 + p16 + p17 + p19
na =~ p2 + p4 + p6 + p7 + p8 + p11 + p13 + p15 + p18 + p20Code
m1 <- cfa(mm, data = d, estimator = "MLR", cluster = "family")
idx <- c("chisq.scaled", "df", "cfi.robust", "rmsea.robust", "srmr")
round(fitMeasures(m1, idx), 3)chisq.scaled df cfi.robust rmsea.robust srmr
1297.385 458.000 0.842 0.059 0.064 Code
round(lavInspect(m1, "cor.lv"), 2) per con pa na
per 1.00
con 0.36 1.00
pa 0.65 0.23 1.00
na -0.26 -0.33 -0.17 1.00Code
std <- standardizedSolution(m1)
std <- std[std$op == "=~", ]
round(tapply(std$est.std, std$lhs, range)[c("per", "con", "pa", "na")] |> do.call(what = rbind), 2) [,1] [,2]
per 0.42 0.76
con 0.16 0.68
pa 0.43 0.73
na 0.49 0.75Code
std[std$est.std < 0.4, c("lhs", "rhs", "est.std")]| lhs | rhs | est.std | |
|---|---|---|---|
| 8 | con | g3 | 0.346920 |
| 12 | con | g11 | 0.156222 |
The fit is modest: CFI 0.842, short of the 0.95 cutoff that Hu and Bentler (1999) suggested, while RMSEA (0.059) and SRMR (0.064) are inside their cutoffs. The other loadings run from 0.42 to 0.76; consistency items 3 and 11 load 0.35 and 0.16.
Code
# Where does the measurement model strain? The largest modification indices: the
# drop in chi-square if one fixed parameter were freed.
mi <- modindices(m1, sort. = TRUE)
head(mi[, c("lhs", "op", "rhs", "mi")], 6)| lhs | op | rhs | mi | |
|---|---|---|---|---|
| 687 | p7 | p20 | 155.21223 | |
| 116 | per | =~ | g11 | 84.60861 |
| 688 | p8 | p11 | 49.90510 | |
| 124 | per | =~ | p16 | 49.62030 |
| 115 | per | =~ | g8 | 47.35594 |
| 551 | p5 | p10 | 44.45408 |
Where is the misfit? The largest modification index is a shared error between PANAS items 7 and 20, scared and afraid (Watson et al., 1988); items 8 and 11 (hostile, irritable) and 5 and 10 (strong, proud) follow. Next is item 11 loading on perseverance, negatively. The strain is among the items: adjectives that share meaning beyond their factor, and a grit item that relates to the other facet. None of it concerns the relations between the constructs, and I read the paths below with that in mind.
Step 2: the paths
Code
# Step 2: replace the four factor correlations with regressions. "~" reads "is
# regressed on"; lavaan correlates the two outcomes' residuals by default.
sm <- paste(mm, "
pa ~ per + con
na ~ per + con")
s1 <- sem(sm, data = d, estimator = "MLR", cluster = "family")
round(rbind(measurement = fitMeasures(m1, idx), structural = fitMeasures(s1, idx)), 3) chisq.scaled df cfi.robust rmsea.robust srmr
measurement 1297.385 458 0.842 0.059 0.064
structural 1297.385 458 0.842 0.059 0.064Code
paths <- standardizedSolution(s1)
paths <- paths[paths$op == "~", c("lhs", "rhs", "est.std", "se", "pvalue")]
paths$est.std <- round(paths$est.std, 2); paths$se <- round(paths$se, 3); paths$pvalue <- round(paths$pvalue, 3)
paths| lhs | rhs | est.std | se | pvalue | |
|---|---|---|---|---|---|
| 33 | pa | per | 0.66 | 0.064 | 0.000 |
| 34 | pa | con | -0.01 | 0.065 | 0.842 |
| 35 | na | per | -0.16 | 0.076 | 0.032 |
| 36 | na | con | -0.27 | 0.072 | 0.000 |
The fit is identical to the measurement model’s, to the last digit: the structural part is saturated, as the builder widget said. Perseverance predicts positive affect strongly (0.66) and lower negative affect a little (−0.16); consistency predicts lower negative affect (−0.27) and nothing about positive affect (−0.01). In words: respondents who say they finish what they begin report feeling more interested, enthusiastic and alert, and those whose interests stay put report somewhat less distress.
What the sum scores would have said
Code
# The same regressions on mean scores, all standardized so the coefficients are on
# the same footing as the standardized latent paths.
alpha <- function(x) { k <- ncol(x); k / (k - 1) * (1 - sum(apply(x, 2, var)) / var(rowSums(x))) }
round(c(per = alpha(d[per]), con = alpha(d[con]), pa = alpha(d[pa]), na = alpha(d[na])), 2) per con pa na
0.74 0.69 0.85 0.87 Code
sc <- as.data.frame(scale(sapply(list(per = per, con = con, pa = pa, na = na),
function(v) rowMeans(d[v]))))
ols <- rbind(coef(lm(pa ~ per + con, sc))[-1], coef(lm(na ~ per + con, sc))[-1])
cbind(paths[, c("lhs", "rhs")], sum_scores = round(as.vector(t(ols)), 2),
latent = paths$est.std)| lhs | rhs | sum_scores | latent | |
|---|---|---|---|---|
| 33 | pa | per | 0.51 | 0.66 |
| 34 | pa | con | 0.03 | -0.01 |
| 35 | na | per | -0.18 | -0.16 |
| 36 | na | con | -0.24 | -0.27 |
With alphas from 0.69 to 0.87, the sum-score regression puts perseverance’s path to positive affect at 0.51, against 0.66 between the factors. The small paths move less, and not all toward zero: perseverance’s path to negative affect is −0.18 on sum scores and −0.16 between factors, the leak the first widget showed. Set the widget to these values (paths 0.66 and 0, correlation 0.36, reliabilities 0.74, 0.69 and 0.85) and it gives 0.51 and 0.04, against the 0.51 and 0.03 printed here.
Twins, and ordered categories
Code
# What the family clustering does to the standard errors of the paths (unstandardized).
s0 <- sem(sm, data = d, estimator = "MLR")
pe1 <- parameterEstimates(s1); pe0 <- parameterEstimates(s0)
data.frame(pe1[pe1$op == "~", c("lhs", "rhs")], est = round(pe1$est[pe1$op == "~"], 2),
se_ignoring_families = round(pe0$se[pe0$op == "~"], 3),
se_clustered = round(pe1$se[pe1$op == "~"], 3))| lhs | rhs | est | se_ignoring_families | se_clustered | |
|---|---|---|---|---|---|
| 33 | pa | per | 0.72 | 0.109 | 0.113 |
| 34 | pa | con | -0.01 | 0.052 | 0.053 |
| 35 | na | per | -0.25 | 0.115 | 0.119 |
| 36 | na | con | -0.30 | 0.086 | 0.087 |
Ignoring the families understates the standard errors, but only slightly here: 0.109 against 0.113 for the largest path (problem 3).
Code
# A check: the items as ordered categories (WLSMV, polychoric correlations).
# Negative-affect items pile up at 1, which is where the two routes could part.
table(negative_affect = unlist(d[na]))negative_affect
1 2 3 4 5
3114 1732 763 382 189 Code
s_ord <- sem(sm, data = d, ordered = TRUE)
po <- standardizedSolution(s_ord)
round(setNames(po$est.std[po$op == "~"], paste(po$lhs, "~", po$rhs)[po$op == "~"]), 2)pa ~ per pa ~ con na ~ per na ~ con
0.67 -0.03 -0.15 -0.27 The negative-affect items are skewed (3,114 of the 6,180 responses are 1), so we refit with the items as ordered categories. The paths barely move (0.67, −0.03, −0.15, −0.27).
Children and adolescents
Do the paths hold for children aged 9–13 and adolescents aged 14–18? A multigroup SEM fits the model in both groups at once and adds constraints: first on the loadings, then on the paths.
Code
# Does the structural model hold for children and for adolescents? Multigroup SEM:
# configural (same model in both), metric (loadings equal), then the four paths equal.
d$ages <- ifelse(d$age < 14, "9-13", "14-18")
table(d$ages, useNA = "ifany") # one respondent has no age and sits this out
14-18 9-13 <NA>
334 283 1 Code
dg <- d[!is.na(d$ages), ]
g_cfg <- sem(sm, data = dg, group = "ages", estimator = "MLR", cluster = "family")
g_met <- sem(sm, data = dg, group = "ages", estimator = "MLR", cluster = "family",
group.equal = "loadings")
g_pth <- sem(sm, data = dg, group = "ages", estimator = "MLR", cluster = "family",
group.equal = c("loadings", "regressions"))
round(rbind(configural = fitMeasures(g_cfg, idx), metric = fitMeasures(g_met, idx),
equal_paths = fitMeasures(g_pth, idx)), 3) chisq.scaled df cfi.robust rmsea.robust srmr
configural 1856.982 916 0.833 0.062 0.073
metric 1920.782 944 0.825 0.062 0.078
equal_paths 1920.164 948 0.826 0.062 0.079Code
lavTestLRT(g_cfg, g_met, g_pth)| Df | AIC | BIC | Chisq | Chisq diff | RMSEA | Df diff | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| g_cfg | 916 | 53031.96 | 53934.64 | 2110.282 | NA | NA | NA | NA |
| g_met | 944 | 53058.53 | 53837.31 | 2192.848 | 62.881554 | 0.0728168 | 28 | 0.0001723 |
| g_pth | 948 | 53052.03 | 53813.11 | 2194.348 | 1.060536 | 0.0000000 | 4 | 0.9004820 |
Code
# Which equal loading strains most? One score test per equality constraint.
st <- lavTestScore(g_met)$uni
pt <- parTable(g_met)
st$constraint <- sapply(st$lhs, function(l) with(pt[pt$plabel == l, ], paste(lhs, op, rhs)))
st <- st[order(-st$X2), c("constraint", "X2", "df", "p.value")]
head(st, 3)| constraint | X2 | df | p.value | |
|---|---|---|---|---|
| 10 | con =~ g11 | 10.761952 | 1 | 0.0010361 |
| 17 | pa =~ p16 | 8.831104 | 1 | 0.0029614 |
| 28 | na =~ p20 | 6.461775 | 1 | 0.0110219 |
Code
# The paths, by age group (unstandardized, under equal loadings).
pg <- parameterEstimates(g_met)
pg <- pg[pg$op == "~", ]
data.frame(path = paste(pg$lhs, "~", pg$rhs), ages = lavInspect(g_met, "group.label")[pg$group],
est = round(pg$est, 2), se = round(pg$se, 2))| path | ages | est | se |
|---|---|---|---|
| pa ~ per | 14-18 | 0.70 | 0.12 |
| pa ~ con | 14-18 | -0.02 | 0.07 |
| na ~ per | 14-18 | -0.25 | 0.15 |
| na ~ con | 14-18 | -0.38 | 0.12 |
| pa ~ per | 9-13 | 0.72 | 0.13 |
| pa ~ con | 9-13 | -0.02 | 0.08 |
| na ~ per | 9-13 | -0.30 | 0.15 |
| na ~ con | 9-13 | -0.23 | 0.11 |
Equal loadings cost something by the test (a scaled Δχ² of 62.9 on 28 df, \(p < 0.001\)), though CFI moves only from 0.833 to 0.825, within the change of 0.01 that Chen (2007) suggests tolerating. The loading under most strain is, again, item 11’s. Given equal loadings, equal paths cost nothing (1.1 on 4 df, \(p = 0.90\)): perseverance predicts positive affect about as strongly at 14–18 as at 9–13 (0.70 and 0.72, unstandardized). That comparison rests on metric invariance, which holds approximately here, not exactly.
Problems
- Attenuation of a slope. Let \(X = T_X + E_X\) and \(Y = \beta T_X + \zeta + E_Y\), with each error uncorrelated with everything else. Show that the population slope of \(Y\) on \(X\) is \(\beta \rho_{XX'}\) and that \(E_Y\) doesn’t bias it. Then find the standardized slope and compare it with Spearman’s formula.
- Factor scores, with a twist. Refit the sum-score regressions with factor scores from the measurement model (
lavPredict(m1)) in place of the sum scores. Are the paths closer to the latent ones or to the sum-score ones? Why doesn’t using factor scores remove the bias, and which kind of factor score for the predictors and which for the outcome would (Skrondal & Laake, 2001)? - Twins. The participants come in pairs. What does clustering do to standard errors, and why did it do so little here? Compute the correlation of twins’ sum scores within families. What would a model that used the pairing as information, rather than a nuisance, look like?
- Design. Draw an SEM for a question you care about, with at least two latent predictors and one latent outcome. Count its structural degrees of freedom. Which path would you most like to read causally, and what design would let you?
- Mend the measurement model. Free the error covariance between PANAS items 7 and 20 (and, if you like, 8 and 11). How much does fit improve, and do the structural paths change? Then drop items 3 and 11 from consistency of interest. Which change matters more for the paths, and which would you defend in a paper?
- Challenge (open). Write down a model equivalent to ours in which affect predicts grit, and one in which a third factor causes both. What evidence (waves, twins, an intervention) could separate them? The PANAS was asked in all three waves, grit in one; I don’t know how far these data can go.
Going further
- Bollen, K. A. (1989). Structural equations with latent variables. Wiley. https://doi.org/10.1002/9781118619179. The standard reference: the model, identification, estimation and fit.
- Jöreskog, K. G. (1970). A general method for estimating a linear structural equation system. ETS Research Bulletin Series, 1970(2). https://doi.org/10.1002/j.2333-8504.1970.tb00783.x. Where factors and paths were first put together and fitted by maximum likelihood.
- 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; its tutorial covers everything used here.
- Anderson, J. C., & Gerbing, D. W. (1988). Structural equation modeling in practice: A review and recommended two-step approach. Psychological Bulletin, 103(3), 411–423. https://doi.org/10.1037/0033-2909.103.3.411. The case for fitting the measurement model first.
- Cole, D. A., & Preacher, K. J. (2014). Manifest variable path analysis: Potentially serious and misleading consequences due to uncorrected measurement error. Psychological Methods, 19(2), 300–315. https://doi.org/10.1037/a0033805. What error does to path models on observed scores, in larger models than ours.
- Westfall, J., & Yarkoni, T. (2016). Statistically controlling for confounding constructs is harder than you think. PLOS ONE, 11(3), e0152719. https://doi.org/10.1371/journal.pone.0152719. The leak in the first widget, and the false positives it produces.
- MacCallum, R. C., Wegener, D. T., Uchino, B. N., & Fabrigar, L. R. (1993). The problem of equivalent models in applications of covariance structure analysis. Psychological Bulletin, 114(1), 185–199. https://doi.org/10.1037/0033-2909.114.1.185. Rules for generating equivalent models, and how many published models have them.
- Bollen, K. A., & Pearl, J. (2013). Eight myths about causality and structural equation models. In S. L. Morgan (Ed.), Handbook of causal analysis for social research (pp. 301–328). Springer. https://doi.org/10.1007/978-94-007-6094-3_15. What the arrows in an SEM do and don’t claim.
- Muthén, B., & Satorra, A. (1995). Complex sample data in structural equation modeling. Sociological Methodology, 25, 267–316. https://doi.org/10.2307/271070. Standard errors for clustered samples, the idea behind
cluster =. - Schuberth, F., Henseler, J., & Dijkstra, T. K. (2018). Confirmatory composite analysis. Frontiers in Psychology, 9, 2541. https://doi.org/10.3389/fpsyg.2018.02541. Composites as a model in their own right.
- 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, 61(5), 745–764. https://doi.org/10.1080/00273171.2026.2645212. Prediction as a complement to fit for the measurement model.
- Muthén, B., & Asparouhov, T. (2016). Multidimensional, multilevel, and multi-timepoint item response modeling. In W. J. van der Linden (Ed.), Handbook of item response theory: Vol. 1. Models (pp. 527–539). Chapman and Hall/CRC. https://doi.org/10.1201/9781315374512. The same ideas with item response models as the measurement model, including longitudinal ones.
- Regressions of ability on covariates inside an item response model come up in Explanatory item response models.
- From the IRW: the CFA walkthrough with lavaan, whose measurement models this lesson extends.
Data sources
Data from the Item Response Warehouse (each table name links to its IRW page); references from IRW biblio.
florida_twins_grit: Hart, S. A., Schatschneider, C., & Taylor, J. E. (2021). Wave 1, Wave 2, Wave 3 Child Survey Measures. LDbase. https://doi.org/10.33009/ldbase.1624481381.d9ec. Licence: CC BY 4.0.florida_twins_panas: Hart, S. A., Schatschneider, C., & Taylor, J. E. (2021). Wave 1, Wave 2, Wave 3 Child Survey Measures. LDbase. https://doi.org/10.33009/ldbase.1624481381.d9ec. Licence: CC BY 4.0.
For instructors
- Source: new material (Ben’s suggestion, #64). EDUC 252 did not cover SEM; the IRW’s lavaan CFA vignette was the starting point.
- Timing: about 50 minutes for the core ideas and widgets; the real-data section is a 30-minute lab (each fit takes seconds locally, and the whole script about a minute).
- Placement: an extension after Factor analysis II, which introduces
lavaan, fit indices and the choice of estimator, and Measurement invariance, which introduces multigroup invariance; this lesson adds~,~~, paths and the two-step approach. - Held back: longitudinal SEM (the Florida twins have three waves of PANAS; problem 6 points at it); twin models that use the pairs to separate genetic and environmental variance; mediation and indirect effects; latent interactions.
- Data: grit and PANAS come from one study and count as one data source. Grit keying was checked against the IRW option text: the perseverance items are reversed and the consistency items are not, unlike the
grittable in Factor analysis II. The family variable is inferred from the id pattern (ids in pairs ending 00 and 01), not from a codebook. The PANAS items are single adjectives from a published scale (Watson et al., 1988); the lesson names six of them. The grit items are described, not reproduced. - Solutions: held back in
solutions/sem.qmd(not published).