library(haven)
library(lavaan)
library(psych)
library(texreg)
library(tidySEM)Reporting lavaan results
APA write-ups, exporting tables with texreg, exporting path diagrams
Learning objectives
By the end of this tutorial you will be able to:
- write an APA-style results paragraph for a CFA, path analysis, SEM, or mediation model
- export lavaan parameter estimates and fit indices to publication-quality tables with
texreg - export path diagrams with
tidySEMas PNG or PDF files - compare nested models with
lavTestLRT()and report the comparison - assemble a complete results section for a dissertation
Setup
Setup code (you have seen this before)
The canonical workflow
By now you have fitted CFA, path, SEM, and mediation models across Tutorials 06–12. This tutorial pulls together the reporting side: how to go from a fitted lavaan object to a polished results section for a dissertation or paper.
The canonical workflow is:
- Fit the model and check convergence
- Extract fit indices with
fitMeasures() - Extract parameter estimates with
summary()orparameterEstimates() - Export a publication-quality table with
texreg - Export a path diagram with
tidySEM - Write the APA results paragraph
We will demonstrate each step using the two-factor CFA from Tutorial 07 and the latent SEM mediation from Tutorial 12.
Refitting the models
First, we refit the two models we will use as running examples:
# Two-factor CFA (Tutorial 07)
model_cfa <- '
PA =~ pa1 + pa2 + pa3
neg =~ na1 + na2 + na3
PA ~~ neg
'
fit_cfa <- cfa(model_cfa, data = ps, estimator = "MLR")
# Latent SEM mediation (Tutorial 12)
model_sem <- '
con =~ p1 + p2 + p3
pa =~ pa1 + pa2 + pa3
creat =~ creat1 + creat2 + creat3
creat ~ b * pa + c * con
pa ~ a * con
indirect := a * b
direct := c
total := c + (a * b)
'
set.seed(1234)
fit_sem <- sem(model_sem, data = ps,
estimator = "ML", se = "boot", bootstrap = 1000)APA-style write-up templates
CFA results template
Confirmatory factor analysis was used to test the hypothesised factor structure of [construct], with [factor 1] measured by [items] and [factor 2] measured by [items]. The model was estimated with robust maximum likelihood (MLR). The model fit the data well, \(\chi^2(df) = value\), \(p = value\), CFI = value, TLI = value, RMSEA = value, 90% CI [lower, upper], SRMR = value. All factor loadings were significant (\(p < .001\)) and substantial (standardised \(\lambda\) range: [low] to [high]). The correlation between [factor 1] and [factor 2] was \(r = value\), \(p < .001\). Cronbach’s alpha was [alpha1] for [factor 1] and [alpha2] for [factor 2], indicating acceptable internal consistency.
SEM mediation results template
Structural equation modelling was used to test whether [mediator] mediated the relationship between [IV] and [DV]. [IV], [mediator], and [DV] were modelled as latent factors measured by [indicators]. The structural paths were estimated with maximum likelihood and [N] bootstrap replications; bias-corrected and accelerated (BCa) confidence intervals are reported for the indirect effect. The model fit the data well, \(\chi^2(df) = value\), CFI = value, TLI = value, RMSEA = value, 90% CI [lower, upper], SRMR = value. The indirect effect of [IV] on [DV] through [mediator] was \(ab = value\), 95% BCa CI [lower, upper], indicating [significant/non-significant] mediation. The direct effect was \(c' = value\), \(p = value\), suggesting [full/partial] mediation. The total effect was \(c = value\), \(p = value\).
Filling in the values
Let us extract the values for the CFA template:
fitMeasures(fit_cfa, c("chisq", "df", "pvalue",
"cfi", "tli",
"rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
"srmr"))## chisq df pvalue cfi tli
## 9.870 8.000 0.274 0.998 0.996
## rmsea rmsea.ci.lower rmsea.ci.upper srmr
## 0.027 0.000 0.073 0.015pe_cfa <- parameterEstimates(fit_cfa, standardized = TRUE)
# Standardised loadings
pe_cfa[pe_cfa$op == "=~", c("rhs", "std.all")]## rhs std.all
## 1 pa1 0.829
## 2 pa2 0.764
## 3 pa3 0.856
## 4 na1 0.848
## 5 na2 0.819
## 6 na3 0.836# Factor correlation
pe_cfa[pe_cfa$op == "~~" & pe_cfa$lhs == "PA" & pe_cfa$rhs == "neg",
c("lhs", "rhs", "std.all", "pvalue")]## lhs rhs std.all pvalue
## 7 PA neg -0.436 0# Cronbach's alpha
psych::alpha(ps[, c("pa1", "pa2", "pa3")])$total$raw_alpha## [1] 0.8482971psych::alpha(ps[, c("na1", "na2", "na3")])$total$raw_alpha## [1] 0.8649448And for the SEM mediation template:
fitMeasures(fit_sem, c("chisq", "df", "pvalue",
"cfi", "tli",
"rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
"srmr"))## chisq df pvalue cfi tli
## 25.528 24.000 0.378 0.999 0.998
## rmsea rmsea.ci.lower rmsea.ci.upper srmr
## 0.014 0.000 0.047 0.028pe_sem <- parameterEstimates(fit_sem, boot.ci.type = "bca.simple")
pe_sem[pe_sem$lhs %in% c("indirect", "direct", "total"),
c("lhs", "est", "ci.lower", "ci.upper", "pvalue")]## lhs est ci.lower ci.upper pvalue
## 25 indirect 0.047 0.011 0.098 0.037
## 26 direct -0.102 -0.241 0.041 0.149
## 27 total -0.055 -0.208 0.091 0.463Exporting tables with texreg
The texreg package exports lavaan results to clean, publication-quality tables in three formats: screen (for checking), HTML (for Word documents), and LaTeX (for PDF documents).
Screen output
For a quick check, use screenreg():
texreg::screenreg(fit_cfa,
summaries = c("N", "AIC", "BIC"),
labels = list(),
single.row = TRUE)##
## ===================================
## Model 1
## -----------------------------------
## PA =~ pa1 1.00 (0.00)
## PA =~ pa2 1.00 (0.06) ***
## PA =~ pa3 1.09 (0.07) ***
## neg =~ na1 1.00 (0.00)
## neg =~ na2 0.94 (0.06) ***
## neg =~ na3 0.98 (0.06) ***
## PA ~~ neg -0.40 (0.07) ***
## pa1 ~~ pa1 0.43 (0.05) ***
## pa2 ~~ pa2 0.68 (0.07) ***
## pa3 ~~ pa3 0.41 (0.06) ***
## na1 ~~ na1 0.35 (0.05) ***
## na2 ~~ na2 0.39 (0.05) ***
## na3 ~~ na3 0.38 (0.04) ***
## PA ~~ PA 0.95 (0.10) ***
## neg ~~ neg 0.90 (0.08) ***
## -----------------------------------
## agfi 0.97
## AIC 5246.25
## BIC 5295.64
## cfi 1.00
## chisq 9.87
## npar 13.00
## rmsea 0.03
## rmsea.conf.high 0.07
## srmr 0.01
## tli 1.00
## converged 1
## ngroups 1
## nobs 330
## norig 343
## nexcluded 13
## ===================================
## *** p < 0.001; ** p < 0.01; * p < 0.05This prints a compact table to the R console with each parameter as a row, its estimate, standard error, and significance stars. The summaries argument controls which fit statistics appear at the bottom.
HTML output (for Word documents)
To insert a lavaan table into a Word document, export it as HTML and copy-paste, or insert the HTML file into Word:
texreg::htmlreg(fit_cfa, file = "tables/table_cfa.html",
summaries = c("N", "AIC", "BIC"),
caption = "Confirmatory Factor Analysis of Positive and Negative Affect",
caption.above = TRUE)Open the HTML file in a web browser, select all, copy, and paste into Word — the table formatting is preserved.
LaTeX output (for PDF documents)
For documents typeset in LaTeX (e.g. via Quarto with a PDF output):
texreg::texreg(fit_cfa, file = "tables/table_cfa.tex",
summaries = c("N", "AIC", "BIC"),
caption = "Confirmatory Factor Analysis of Positive and Negative Affect",
caption.above = TRUE,
label = "tab:cfa")Include the table in a LaTeX or Quarto document with \input{tables/table_cfa.tex} or via Quarto’s $include$ mechanism.
Comparing multiple models
texreg can display several models side by side, which is useful for comparing nested models (e.g. a first-order CFA vs. a second-order CFA, or a path model vs. a full SEM):
texreg::screenreg(list(fit_cfa, fit_cfa),
custom.model.names = c("First-order", "Second-order"),
summaries = c("N", "AIC", "BIC"))
Note
texreg and lavaan
texreg works with lavaan objects out of the box — it knows how to extract parameter estimates, fit indices, and sample size. For mediation models with := defined parameters, the defined parameters (indirect, direct, total) appear in the table alongside the regression coefficients.
For more complex table layouts, consider semTable (from the semTable package) as an alternative — it is specifically designed for SEM and produces slightly different formatting. We recommend texreg as the default because it handles both lm and lavaan objects, making it useful across your entire dissertation.
Exporting path diagrams
With tidySEM
To save a path diagram as a PNG or PDF, use tidySEM’s graph functions with a graphics device:
# Save as PNG
png("figures/sem_mediation.png", width = 1200, height = 800, res = 150)
tidySEM::graph_sem(fit_sem, std = TRUE)
dev.off()
# Save as PDF (vector graphics, publication quality)
pdf("figures/sem_mediation.pdf", width = 8, height = 6)
tidySEM::graph_sem(fit_sem, std = TRUE)
dev.off()
TipTry it
Fit the two-factor CFA from Tutorial 07 and export its path diagram as a PNG using the code above. Open the PNG file and check that the loadings and factor correlation are visible on the arrows.
CautionSolution
fit_cfa <- cfa(model_cfa, data = ps, estimator = "MLR")
png("figures/cfa_pa_na.png", width = 1200, height = 800, res = 150)
tidySEM::graph_sem(fit_cfa, std = TRUE)
dev.off()The PNG file in figures/cfa_pa_na.png should show the two latent factors (PA, neg) with their indicators and the standardised loadings on the arrows, plus the factor correlation on the double-headed arrow between them.
Comparing nested models
lavTestLRT() for MLR models
When comparing nested models fitted with MLR, use lavTestLRT() rather than anova() — the Satorra–Bentler scaled chi-square difference test is the correct test for robust estimators:
# First-order CFA (Tutorial 07)
fit_1st <- cfa(model_cfa, data = ps, estimator = "MLR")
# Second-order CFA (Tutorial 08)
model_2nd <- '
PA =~ pa1 + pa2 + pa3
neg =~ na1 + na2 + na3
affect =~ PA + neg
'
fit_2nd <- cfa(model_2nd, data = ps, estimator = "MLR")
lavTestLRT(fit_1st, fit_2nd)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_2nd 7 5248.3 5301.4 9.87
## fit_1st 8 5246.3 5295.6 9.87 0 1The output reports the scaled chi-square difference, its degrees of freedom, and the p-value. A non-significant result means the more restrictive model (here, the second-order) fits as well as the less restrictive model — the restriction is justified.
anova() for ML models
For models fitted with plain ML (e.g. bootstrap models), anova() is appropriate:
anova(fit_1, fit_2) # both fitted with estimator = "ML"
WarningUse the right test
MLRmodels →lavTestLRT()(scaled difference test)MLmodels →anova()(standard likelihood-ratio test)- Bootstrap models → no standard difference test; compare bootstrap CIs instead
Using the wrong test gives incorrect p-values. When in doubt, lavTestLRT() is the safer choice because it works for both ML and MLR.
A complete results section
Putting it all together — here is a complete results section for the latent SEM mediation model from Tutorial 12, using the values extracted above:
Model description
Structural equation modelling was used to test whether positive affect (PA) mediated the relationship between job control and creativity. Job control was modelled as a latent factor measured by three parcels (p1, p2, p3) constructed from the six job-control items. PA and creativity were modelled as latent factors measured by their three respective items. The structural paths were estimated with maximum likelihood and 1,000 bootstrap replications; bias-corrected and accelerated (BCa) confidence intervals are reported for the indirect effect.
Model fit
The model fit the data well (see Table 1), \(\chi^2(24) = ...\), CFI = …, TLI = …, RMSEA = …, 90% CI […, …], SRMR = …. All factor loadings were significant (\(p < .001\)) and substantial, confirming that the measurement model was adequate.
Mediation results
The indirect effect of job control on creativity through PA was \(ab = ...\), 95% BCa CI […, …]. Because the confidence interval did not include zero, the indirect effect was significant, indicating that PA mediated the relationship between job control and creativity. The direct effect of job control on creativity was \(c' = ...\), \(p = ...\), suggesting [full/partial] mediation. The total effect was \(c = ...\), \(p = ...\).
Path diagram
[Insert path diagram here — exported with tidySEM::graph_sem() as shown above]
Table
[Insert table here — exported with texreg::htmlreg() as shown above]
Common errors
WarningReporting a just-identified model’s fit
Do not report “excellent fit” for a just-identified model (df = 0). The fit is guaranteed by the model’s structure. Report the parameter estimates and \(R^2\) values, and note that the model is just-identified.
Using anova() with MLR models
anova() does not apply the Satorra–Bentler scaling correction. Use lavTestLRT() for models fitted with MLR. The two functions give different p-values; only lavTestLRT() is correct for robust estimators.
Forgetting to set the bootstrap seed
If you re-run a bootstrap model without set.seed(), the CIs will differ. Always set the seed, and report the seed in your methods section for full reproducibility.
Exercises
Fit the two-factor CFA from Tutorial 07. Extract the fit indices and the standardised loadings with
parameterEstimates(). Write a complete APA-style results paragraph using the CFA template above, filling in every value.Fit the latent SEM mediation model from Tutorial 12. Export the parameter estimates to an HTML table with
texreg::htmlreg(). Save the path diagram as a PNG withtidySEM::graph_sem().Fit a first-order two-factor CFA and a second-order CFA (both on PA and neg, as in Tutorial 08). Compare them with
lavTestLRT(). Write a one-paragraph summary of the comparison, reporting the scaled chi-square difference, its degrees of freedom, and its p-value.Fit the path mediation model from Tutorial 12 (manifest variables,
MLwithse = "boot"). Extract the BCa CI for the indirect effect. Write an APA-style results paragraph for this model, following the SEM mediation template.
CautionSolutions
Exercise 1
fit_cfa <- cfa(model_cfa, data = ps, estimator = "MLR")
fitMeasures(fit_cfa, c("chisq", "df", "pvalue",
"cfi", "tli",
"rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "srmr"))## chisq df pvalue cfi tli
## 9.870 8.000 0.274 0.998 0.996
## rmsea rmsea.ci.lower rmsea.ci.upper srmr
## 0.027 0.000 0.073 0.015pe <- parameterEstimates(fit_cfa, standardized = TRUE)
pe[pe$op == "=~", c("rhs", "std.all")]## rhs std.all
## 1 pa1 0.829
## 2 pa2 0.764
## 3 pa3 0.856
## 4 na1 0.848
## 5 na2 0.819
## 6 na3 0.836psych::alpha(ps[, c("pa1", "pa2", "pa3")])$total$raw_alpha## [1] 0.8482971psych::alpha(ps[, c("na1", "na2", "na3")])$total$raw_alpha## [1] 0.8649448Fill in each value in the CFA template. Example: “\(\chi^2(8) = 8.79\), \(p = .360\), CFI = 1.00, TLI = 1.00, RMSEA = 0.02, 90% CI [0.00, 0.06], SRMR = 0.01.”
Exercise 2
texreg::htmlreg(fit_sem, file = "tables/table_sem.html",
caption = "Latent SEM Mediation Model",
caption.above = TRUE)
png("figures/sem_mediation.png", width = 1200, height = 800, res = 150)
tidySEM::graph_sem(fit_sem, std = TRUE)
dev.off()Exercise 3
fit_1st <- cfa(model_cfa, data = ps, estimator = "MLR")
fit_2nd <- cfa(model_2nd, data = ps, estimator = "MLR")
lavTestLRT(fit_1st, fit_2nd)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_2nd 7 5248.3 5301.4 9.87
## fit_1st 8 5246.3 5295.6 9.87 0 1Write: “The second-order model was compared to the first-order model using the Satorra–Bentler scaled chi-square difference test. The second-order restriction did not significantly worsen fit, \(\Delta\chi^2(...) = ...\), \(p = ...\), indicating that the more parsimonious second-order structure was supported.”
Exercise 4
set.seed(1234)
ps$con <- (ps$con1 + ps$con3 + ps$con6 + ps$con8 + ps$con9 + ps$con10) / 6
ps$pa <- (ps$pa1 + ps$pa2 + ps$pa3) / 3
ps$creat <- (ps$creat1 + ps$creat2 + ps$creat3) / 3
model_path <- '
creat ~ b * pa + c * con
pa ~ a * con
indirect := a * b
direct := c
total := c + (a * b)
'
fit_path <- sem(model_path, data = ps,
estimator = "ML", se = "boot", bootstrap = 1000)
pe_path <- parameterEstimates(fit_path, boot.ci.type = "bca.simple")
pe_path[pe_path$lhs %in% c("indirect", "direct", "total"),
c("lhs", "est", "ci.lower", "ci.upper", "pvalue")]## lhs est ci.lower ci.upper pvalue
## 7 indirect 0.042 0.012 0.085 0.022
## 8 direct -0.091 -0.226 0.043 0.186
## 9 total -0.049 -0.195 0.096 0.504Fill in the SEM mediation template with the extracted values.
Summary & next tutorial
You can now write APA-style results paragraphs for CFA, path analysis, SEM, and mediation models; export parameter estimates and fit indices to publication-quality tables with texreg; export path diagrams with tidySEM; and compare nested models with the correct statistical test. This completes the core SEM workflow — you have everything you need for a dissertation using lavaan. In the remaining tutorials we cover advanced topics: multilevel models, growth curve models, and moderation, which extend the lavaan toolkit to nested, longitudinal, and interaction-based analyses.
Further reading
- Rosseel (2012) — the lavaan reference paper.
- Beaujean (2014) — a lavaan-based walkthrough including reporting.
- Kelloway (2015) — the textbook basis for the original teaching materials; includes chapters on reporting SEM results.
- Hu & Bentler (1999) — the source of the fit-index cutoff criteria used in APA write-ups.
References
Beaujean, A. A. (2014). Latent variable modeling using r: A step-by-step guide. Routledge.
Hu, L., & Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling: A Multidisciplinary Journal, 6(1), 1–55. https://doi.org/10.1080/10705519909540118
Kelloway, E. K. (2015). Using mplus for structural equation modeling (2nd ed.). SAGE Publications.
Rosseel, Y. (2012). Lavaan: An r package for structural equation modeling. Journal of Statistical Software, 48(2), 1–36. https://doi.org/10.18637/jss.v048.i02