CFA — practical

Two-factor CFA on NA/PA; fit, parameter estimates, modification indices

Author
Affiliation

George Michaelides

Norwich Business School, University of East Anglia

Published

July 7, 2026

Learning objectives

By the end of this tutorial you will be able to:

  • specify and fit a two-factor CFA in lavaan using the PS dataset
  • read every section of the summary(fit, standardized = TRUE, fit.measures = TRUE) output
  • extract the canonical fit indices with fitMeasures() and judge them against the Hu and Bentler cutoffs
  • compute Cronbach’s alpha for a scale with psych::alpha()
  • request and interpret modification indices with modindices()
  • report CFA results in APA format

Setup

Setup code (you have seen this before) 
library(haven)
library(lavaan)
library(psych)
library(tidySEM)

From Tutorial 06 onward, the setup block — loading packages, importing the PS dataset, and converting missing values — is collapsed behind a <details> tag. Click to expand it if you need to copy the code. The block is the same in every tutorial that uses the PS data.

The model

We fit the two-factor CFA on positive and negative affect that we introduced in Tutorial 06:

  • Positive affect (PA) measured by pa1, pa2, pa3
  • Negative affect (neg) measured by na1, na2, na3
  • The two factors are allowed to correlate

The lavaan syntax:

model_pa_na <- '
  PA  =~ pa1 + pa2 + pa3
  neg =~ na1 + na2 + na3
  PA  ~~ neg
'

We use the marker-variable method (lavaan’s default) to set each factor’s scale: the first loading of each factor is fixed to 1. We estimate with MLR — our default throughout the series.

Fitting the model

fit_pa_na <- cfa(model_pa_na, data = ps, estimator = "MLR", missing = "ml")

The model converges silently — no warnings, which means lavaan found a stable solution. If it had not converged, we would see a warning and would need to diagnose the problem (see the Common errors section at the end).

Reading the full output

We call summary() with both standardized = TRUE and fit.measures = TRUE to get the complete picture: fit indices at the top, parameter estimates below, with both unstandardised and standardised values.

summary(fit_pa_na, standardized = TRUE, fit.measures = TRUE)
## lavaan 0.6-21 ended normally after 32 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        19
## 
##                                                   Used       Total
##   Number of observations                           342         343
##   Number of missing patterns                         6            
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                10.298       8.789
##   Degrees of freedom                                 8           8
##   P-value (Chi-square)                           0.245       0.360
##   Scaling correction factor                                  1.172
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Model Test Baseline Model:
## 
##   Test statistic                              1005.424     803.338
##   Degrees of freedom                                15          15
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.252
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.998       0.999
##   Tucker-Lewis Index (TLI)                       0.996       0.998
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            0.999
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -2706.285   -2706.285
##   Scaling correction factor                                  1.045
##       for the MLR correction                                      
##   Loglikelihood unrestricted model (H1)      -2701.136   -2701.136
##   Scaling correction factor                                  1.082
##       for the MLR correction                                      
##                                                                   
##   Akaike (AIC)                                5450.570    5450.570
##   Bayesian (BIC)                              5523.431    5523.431
##   Sample-size adjusted Bayesian (SABIC)       5463.159    5463.159
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.029       0.017
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.074       0.064
##   P-value H_0: RMSEA <= 0.050                    0.733       0.846
##   P-value H_0: RMSEA >= 0.080                    0.026       0.008
##                                                                   
##   Robust RMSEA                                               0.014
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.073
##   P-value H_0: Robust RMSEA <= 0.050                         0.784
##   P-value H_0: Robust RMSEA >= 0.080                         0.028
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.013       0.013
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   PA =~                                                                 
##     pa1               1.000                               0.976    0.814
##     pa2               1.017    0.069   14.774    0.000    0.992    0.772
##     pa3               1.087    0.076   14.267    0.000    1.061    0.852
##   neg =~                                                                
##     na1               1.000                               0.968    0.846
##     na2               0.956    0.059   16.304    0.000    0.925    0.835
##     na3               0.975    0.056   17.562    0.000    0.944    0.830
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   PA ~~                                                                 
##     neg              -0.433    0.073   -5.923    0.000   -0.459   -0.459
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .pa1               3.632    0.065   56.018    0.000    3.632    3.029
##    .pa2               3.378    0.070   48.492    0.000    3.378    2.630
##    .pa3               3.451    0.068   50.920    0.000    3.451    2.771
##    .na1               2.764    0.062   44.675    0.000    2.764    2.417
##    .na2               2.704    0.060   45.042    0.000    2.704    2.441
##    .na3               2.608    0.062   42.191    0.000    2.608    2.293
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .pa1               0.485    0.079    6.132    0.000    0.485    0.338
##    .pa2               0.666    0.069    9.583    0.000    0.666    0.404
##    .pa3               0.426    0.063    6.746    0.000    0.426    0.275
##    .na1               0.372    0.053    7.062    0.000    0.372    0.284
##    .na2               0.371    0.051    7.333    0.000    0.371    0.303
##    .na3               0.402    0.050    8.079    0.000    0.402    0.311
##     PA                0.952    0.112    8.512    0.000    1.000    1.000
##     neg               0.937    0.090   10.367    0.000    1.000    1.000

Estimator and optimisation

The header tells us the estimator (MLR), the optimisation method (NLMINB), the number of iterations (typically 20–40 for a well-behaved model), and the number of observations used. Because we passed missing = "ml", lavaan uses FIML (full-information maximum likelihood): 342 of 343 participants are retained even if some have missing values on individual items. (One case with no data on any of the six items cannot contribute and is dropped.)

Model fit indices

The Model Fit Information section reports the canonical set:

  • Chi-Square Test of Model Fit: \(\chi^2\) = 8.789 (scaled), df = 8, p = .360. The \(\chi^2\) is non-significant, meaning we cannot reject the null hypothesis that the model-implied covariance matrix equals the observed one. This is good news for our model. The Scaling correction factor is reported because we are using MLR — the \(\chi^2\) has been corrected for non-normality.

  • RMSEA: 0.017, 90% CI [0.000, 0.064]. The point estimate is well below .05, and the upper bound of the CI is below .08. This indicates close fit.

  • CFI/TLI: CFI = 0.999, TLI = 0.998. Both are at the maximum, well above the .95 threshold.

  • SRMR: 0.013. Well below the .08 threshold.

All five indices agree: the model fits the data very well. We will report this in APA format at the end of the tutorial.

NoteWhy is the \(\chi^2\) starred?

The asterisk next to the \(\chi^2\) value in MLR output indicates that the standard \(\chi^2\) difference test does not apply directly to robust estimators. If you need to compare nested models fitted with MLR, use lavTestLRT() (the scaled likelihood-ratio test) rather than anova(). We cover this in Tutorial 13.

Parameter estimates

The Parameter Estimates section has several subsections. Let us walk through each.

Latent variables (loadings)

Latent Variables:
                   Estimate    Std.Err  z value  P(>|z|)   std.all
  PA =~
    pa1              1.000                               0.814
    pa2              1.017       0.069   14.774    0.000    0.772
    pa3              1.087       0.076   14.267    0.000    0.852
  neg =~
    na1              1.000                               0.846
    na2              0.956       0.059   16.304    0.000    0.835
    na3              0.975       0.056   17.562    0.000    0.830

The Estimate column gives the unstandardised factor loading (\(\lambda\)). The first loading of each factor (pa1 and na1) is fixed to 1.000 — the marker-variable method. These fixed loadings have no standard error or p-value, which is why their cells are blank.

The std.all column gives the fully standardised loadings — the correlation between each indicator and its factor. All six loadings are high (0.77–0.85), indicating that each item is a strong indicator of its underlying factor. A loading below .40 would be concerning; below .30 suggests the item is not measuring the factor well.

All freely-estimated loadings have p < .001, meaning each is significantly different from zero.

Covariances

Covariances:
                   Estimate    Std.Err  z value  P(>|z|)   std.all
  PA ~~
    neg             -0.433       0.073   -5.923    0.000   -0.459

The covariance between PA and neg is -0.433. The std.all value (-0.459) is the correlation between the two factors — positive and negative affect are moderately negatively related, which makes theoretical sense. The p < .001 tells us this correlation is significantly different from zero.

Intercepts

Intercepts:
                   Estimate    Std.Err  z value  P(>|z|)   std.all
   .pa1              3.632       0.065   56.018    0.000    3.029
   .pa2              3.378       0.070   48.492    0.000    2.630
     ...

The intercepts are the estimated means of the observed variables. They are not usually of substantive interest in CFA — they tell you the average response on each item — but they are part of the model. The dot prefix (.pa1) indicates that this is an intercept for an observed variable that is dependent on a latent variable.

Variances

Variances:
                   Estimate    Std.Err  z value  P(>|z|)   std.all
     PA               0.952       0.112    8.512    0.000    1.000
     neg              0.937       0.090   10.367    0.000    1.000
    .pa1              0.485       0.079    6.132    0.000    0.338
    .pa2              0.666       0.069    9.583    0.000    0.404
     ...

The first two rows are the latent variances (\(\psi\)) — the variance of each factor. The std.all is 1.000 because a variable’s standardised variance is always 1.

The remaining rows are the residual variances (\(\theta\)) — the portion of each indicator’s variance not explained by the factor. The std.all for a residual variance is the proportion of the indicator’s variance that is not explained by the factor. For pa1, that is 0.338, meaning \(1 - 0.338 = 0.662\) (or 66.2%) of the variance in pa1 is explained by the PA factor. This is the indicator’s \(R^2\).

R-squared

To see the \(R^2\) values explicitly, re-run summary() with rsquare = TRUE:

summary(fit_pa_na, standardized = TRUE, rsquare = TRUE, fit.measures = FALSE)
## lavaan 0.6-21 ended normally after 32 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        19
## 
##                                                   Used       Total
##   Number of observations                           342         343
##   Number of missing patterns                         6            
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                10.298       8.789
##   Degrees of freedom                                 8           8
##   P-value (Chi-square)                           0.245       0.360
##   Scaling correction factor                                  1.172
##     Yuan-Bentler correction (Mplus variant)                       
## 
## Parameter Estimates:
## 
##   Standard errors                             Sandwich
##   Information bread                           Observed
##   Observed information based on                Hessian
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   PA =~                                                                 
##     pa1               1.000                               0.976    0.814
##     pa2               1.017    0.069   14.774    0.000    0.992    0.772
##     pa3               1.087    0.076   14.267    0.000    1.061    0.852
##   neg =~                                                                
##     na1               1.000                               0.968    0.846
##     na2               0.956    0.059   16.304    0.000    0.925    0.835
##     na3               0.975    0.056   17.562    0.000    0.944    0.830
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   PA ~~                                                                 
##     neg              -0.433    0.073   -5.923    0.000   -0.459   -0.459
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .pa1               3.632    0.065   56.018    0.000    3.632    3.029
##    .pa2               3.378    0.070   48.492    0.000    3.378    2.630
##    .pa3               3.451    0.068   50.920    0.000    3.451    2.771
##    .na1               2.764    0.062   44.675    0.000    2.764    2.417
##    .na2               2.704    0.060   45.042    0.000    2.704    2.441
##    .na3               2.608    0.062   42.191    0.000    2.608    2.293
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .pa1               0.485    0.079    6.132    0.000    0.485    0.338
##    .pa2               0.666    0.069    9.583    0.000    0.666    0.404
##    .pa3               0.426    0.063    6.746    0.000    0.426    0.275
##    .na1               0.372    0.053    7.062    0.000    0.372    0.284
##    .na2               0.371    0.051    7.333    0.000    0.371    0.303
##    .na3               0.402    0.050    8.079    0.000    0.402    0.311
##     PA                0.952    0.112    8.512    0.000    1.000    1.000
##     neg               0.937    0.090   10.367    0.000    1.000    1.000
## 
## R-Square:
##                    Estimate
##     pa1               0.662
##     pa2               0.596
##     pa3               0.725
##     na1               0.716
##     na2               0.697
##     na3               0.689

The R-Square section at the bottom lists each indicator’s \(R^2\) — the proportion of its variance explained by its factor. Values above .50 (i.e. standardised loadings above .71) are considered good; all six of our indicators exceed this threshold.

Extracting fit indices cleanly

For reporting, we extract just the canonical set with a named fitMeasures() call:

fit_pa_na_indices <- fitMeasures(fit_pa_na,
                                 c("chisq", "df", "pvalue",
                                   "cfi", "tli",
                                   "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                                   "srmr"))
fit_pa_na_indices
##          chisq             df         pvalue            cfi            tli 
##         10.298          8.000          0.245          0.998          0.996 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##          0.029          0.000          0.074          0.013

This gives a named numeric vector — easy to round and drop into a report. We will use this exact call in every tutorial that fits a model.

Extracting parameter estimates as a data frame

For building tables or exporting to a report, parameterEstimates() returns the estimates as a data frame — one row per parameter:

pe <- parameterEstimates(fit_pa_na, standardized = TRUE)
head(pe[, c("lhs", "op", "rhs", "est", "se", "z", "pvalue", "std.all")], 10)
##    lhs op rhs    est    se      z pvalue std.all
## 1   PA =~ pa1  1.000 0.000     NA     NA   0.814
## 2   PA =~ pa2  1.017 0.069 14.774      0   0.772
## 3   PA =~ pa3  1.087 0.076 14.267      0   0.852
## 4  neg =~ na1  1.000 0.000     NA     NA   0.846
## 5  neg =~ na2  0.956 0.059 16.304      0   0.835
## 6  neg =~ na3  0.975 0.056 17.562      0   0.830
## 7   PA ~~ neg -0.433 0.073 -5.923      0  -0.459
## 8  pa1 ~~ pa1  0.485 0.079  6.132      0   0.338
## 9  pa2 ~~ pa2  0.666 0.069  9.583      0   0.404
## 10 pa3 ~~ pa3  0.426 0.063  6.746      0   0.275

Each row is one parameter. The lhs, op, and rhs columns describe the relationship (e.g. PA =~ pa1), est is the unstandardised estimate, se the standard error, z the Wald statistic, pvalue the two-tailed p-value, and std.all the fully standardised estimate. We will use this data frame in Tutorial 13 when we export results to publication-quality tables.

Cronbach’s alpha

Before or alongside a CFA, it is good practice to report the internal consistency of each scale. Cronbach’s alpha is the standard coefficient. The psych::alpha() function computes it:

psych::alpha(ps[, c("pa1", "pa2", "pa3")])
## 
## Reliability analysis   
## Call: psych::alpha(x = ps[, c("pa1", "pa2", "pa3")])
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
##       0.85      0.85    0.79      0.66 5.7 0.014  3.5 1.1     0.66
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.82  0.85  0.87
## Duhachek  0.82  0.85  0.88
## 
##  Reliability if an item is dropped:
##     raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## pa1      0.79      0.79    0.66      0.66 3.8    0.022    NA  0.66
## pa2      0.81      0.81    0.68      0.68 4.3    0.021    NA  0.68
## pa3      0.76      0.77    0.63      0.63 3.3    0.025    NA  0.63
## 
##  Item statistics 
##       n raw.r std.r r.cor r.drop mean  sd
## pa1 342  0.88  0.88  0.78   0.71  3.6 1.2
## pa2 338  0.87  0.87  0.75   0.70  3.4 1.3
## pa3 332  0.89  0.89  0.81   0.74  3.5 1.2
## 
## Non missing response frequency for each item
##        1    2    3    4    5    6 miss
## pa1 0.04 0.15 0.23 0.34 0.20 0.04 0.00
## pa2 0.07 0.21 0.25 0.26 0.18 0.04 0.01
## pa3 0.05 0.18 0.30 0.26 0.16 0.05 0.03

The output reports raw_alpha (Cronbach’s alpha based on covariances) and std.alpha (standardised alpha). Values above .70 are considered acceptable; above .80, good. The PA scale’s alpha should be comfortably above .70 given the high loadings we observed.

The alpha() output also includes:

  • G6 — Guttman’s lambda 6, an alternative lower bound on reliability
  • average_r — the average inter-item correlation
  • Items — a table showing what would happen to alpha if each item were dropped. If dropping an item increases alpha, that item may be weakening the scale.
TipTry it

Compute Cronbach’s alpha for the negative-affect scale (na1, na2, na3). Is the reliability similar to the PA scale? Would dropping any item improve alpha?

psych::alpha(ps[, c("na1", "na2", "na3")])
## 
## Reliability analysis   
## Call: psych::alpha(x = ps[, c("na1", "na2", "na3")])
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean sd median_r
##       0.86      0.87    0.82      0.69 6.8 0.013  2.7  1     0.69
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.84  0.86  0.89
## Duhachek  0.84  0.86  0.89
## 
##  Reliability if an item is dropped:
##     raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## na1      0.81      0.82    0.69      0.69 4.4    0.021    NA  0.69
## na2      0.81      0.81    0.69      0.69 4.4    0.021    NA  0.69
## na3      0.82      0.83    0.71      0.71 4.8    0.020    NA  0.71
## 
##  Item statistics 
##       n raw.r std.r r.cor r.drop mean  sd
## na1 341  0.90  0.89  0.81   0.75  2.8 1.1
## na2 339  0.89  0.89  0.81   0.75  2.7 1.1
## na3 332  0.89  0.89  0.80   0.73  2.6 1.1
## 
## Non missing response frequency for each item
##        1    2    3    4    5    6 miss
## na1 0.16 0.26 0.31 0.21 0.05 0.01 0.01
## na2 0.14 0.30 0.34 0.16 0.05 0.01 0.01
## na3 0.17 0.32 0.30 0.17 0.04 0.01 0.03

Compare the raw_alpha value to the PA scale’s alpha. The “Item Dropped” column in the Item Statistics table shows what alpha would be if each item were removed. If none of the values exceeds the overall alpha, all three items are contributing positively to the scale’s reliability.

Modification indices

Modification indices (MI) estimate how much the model \(\chi^2\) would improve if a currently-fixed parameter were freely estimated. An MI above 3.84 corresponds to \(\alpha = .05\); above 6.63, to \(\alpha = .01\).

modindices(fit_pa_na, minimum.value = 3.84, sort. = TRUE)
##    lhs op rhs    mi   epc sepc.lv sepc.all sepc.nox
## 36 pa2 ~~ na1 5.840 0.089   0.089    0.178    0.178
## 40 pa3 ~~ na2 4.785 0.071   0.071    0.179    0.179

The output lists every fixed parameter whose MI exceeds 3.84, sorted from largest to smallest. The columns are:

  • lhs, op, rhs — the parameter that would be freed (e.g. pa2 ~~ na1 means “free the residual covariance between pa2 and na1”)
  • mi — the modification index (estimated drop in \(\chi^2\) if this parameter were freed)
  • epc — the expected parameter change (what the parameter’s value would be if freed)
  • sepc.lv, sepc.all, sepc.nox — standardised versions of the EPC
WarningUse modification indices cautiously

Modification indices are data-driven suggestions. They tell you where the sample covariance matrix departs from the model, but they do not tell you why. Freeing a parameter because it has a high MI can improve fit in your sample while making the model worse in the population — this is overfitting. Always have a theoretical reason for freeing a parameter, and report any model changes you make transparently. We recommend using MIs as diagnostic information, not as an automatic model-building tool.

In our two-factor CFA, the modification indices are all small (the largest is below 5), which is consistent with the model’s already-excellent fit. There is no compelling reason to free any additional parameter.

Drawing the path diagram

tidySEM::graph_sem() produces a path diagram directly from the fitted lavaan object:

tidySEM::graph_sem(fit_pa_na, std = TRUE)

Path diagram for the two-factor CFA of positive and negative affect, with standardised loadings.

The diagram shows:

  • Two circles (latent variables: PA, neg)
  • Six squares (observed indicators)
  • Arrows from each factor to its indicators, labelled with the standardised loadings
  • A double-headed arrow between the two factors, labelled with their correlation (-0.46)

This is the same diagrammatic convention we used in Tutorial 06’s concept diagram, but now it is generated from the fitted model — the numbers on the arrows are the actual estimated loadings, not placeholders.

NotetidySEM vs hand-authored diagrams

Throughout this series we use two diagramming tools:

  • tidySEM::graph_sem(fit) — generates a diagram from a fitted lavaan object. The numbers are the actual estimates. Use this when you want to show the results of a fitted model.
  • DiagrammeR::grViz() — hand-authored diagrams for illustrating concepts before a model is fitted. Use this when you want to explain a model structure without tying it to specific estimates.

We used DiagrammeR for the concept diagrams in Tutorials 05 and 06; from Tutorial 07 onward we use tidySEM whenever we want to visualise a fitted model.

Writing up the results

An APA-style write-up of this CFA would look something like:

Confirmatory factor analysis was used to test the hypothesised two-factor structure of affect, with positive affect (PA) measured by three items (pa1, pa2, pa3) and negative affect (neg) measured by three items (na1, na2, na3). The model was estimated with robust maximum likelihood (MLR) to correct for non-normality, using full-information maximum likelihood for missing data. The model fit the data well, \(\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. All factor loadings were significant (\(p < .001\)) and substantial (standardised \(\lambda\) range: .77–.85). The correlation between PA and neg was \(r = -.46\), \(p < .001\). Cronbach’s alpha was .85 for PA and .86 for neg, indicating good internal consistency for both scales.

NoteReporting conventions

In APA format:

  • Report \(\chi^2\) with degrees of freedom and p-value: \(\chi^2(8) = 8.79\), \(p = .360\)
  • Report fit indices to two decimal places: CFI = 1.00, RMSEA = 0.00
  • Report RMSEA with its 90% CI: RMSEA = 0.02, 90% CI [0.00, 0.06]
  • Report standardised loadings with the symbol \(\lambda\) or just “standardised loading”
  • Report correlations with \(r\): \(r = -.46\)
  • Suppress leading zeros for values that cannot exceed 1: \(r = -.46\), not \(r = -0.46\); \(p = .001\), not \(p = 0.001\)

Common errors

Warninglavaan WARNING: model has not converged

The estimation algorithm could not find a stable solution. Possible causes: model misspecification, underidentification, empirical underidentification (two items nearly identical), or a Heywood case. Check summary(fit) for negative variances or loadings above 1, check that every factor has at least three indicators, and check the correlation matrix of your indicators for near-perfect collinearity.

lavaan WARNING: some estimated lv variances are negative

A Heywood case — a negative variance, which is impossible. Usually indicates that the model is misspecified or the sample is too small. Try constraining the problematic variance to a small positive value, or respecify the model.

Very high modification indices on residual covariances

If modindices() suggests large residual covariances among items with similar wording, you may have common method variance. Consider freeing one or two theoretically-justified residual covariances, but report the change transparently.

Exercises

  1. Fit a two-factor CFA with job control (measured by con1, con3, con6, con8, con9, con10) and social support (measured by sup2, sup4, sup5, sup7) as the two factors. Use MLR. How many degrees of freedom does this model have? Does it fit well?

  2. For the job-control factor in Exercise 1, compute Cronbach’s alpha. Would dropping any item improve alpha?

  3. Extract the standardised loadings from Exercise 1 as a data frame using parameterEstimates(fit, standardized = TRUE). Filter the result to show only the factor loadings (rows where op == "=~"). Which item has the highest standardised loading on each factor?

  4. Run modindices() on the model from Exercise 1. Are there any modification indices above 3.84? If so, what parameter would they suggest freeing, and is there a theoretical justification for it?

Exercise 1

model_con_sup <- '
  con =~ con1 + con3 + con6 + con8 + con9 + con10
  sup =~ sup2 + sup4 + sup5 + sup7
  con ~~ sup
'
fit_con_sup <- cfa(model_con_sup, data = ps, estimator = "MLR")
fitMeasures(fit_con_sup, c("chisq", "df", "pvalue",
                           "cfi", "tli",
                           "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                           "srmr"))
##          chisq             df         pvalue            cfi            tli 
##         44.054         34.000          0.116          0.990          0.986 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##          0.029          0.000          0.052          0.034

With \(q = 10\) observed variables, the covariance matrix has \(10 \times 11 / 2 = 55\) elements. Free parameters: 5 loadings on con (first fixed to 1) + 3 loadings on sup (first fixed to 1) + 2 latent variances + 1 covariance + 10 residual variances = 21. So \(\text{df} = 55 - 21 = 34\). Judge the fit indices against the Hu and Bentler cutoffs.

Exercise 2

psych::alpha(ps[, c("con1", "con3", "con6", "con8", "con9", "con10")])
## 
## Reliability analysis   
## Call: psych::alpha(x = ps[, c("con1", "con3", "con6", "con8", "con9", 
##     "con10")])
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.82      0.83    0.81      0.45 4.9 0.015  3.8 0.68     0.45
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.79  0.82  0.85
## Duhachek  0.80  0.82  0.85
## 
##  Reliability if an item is dropped:
##       raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
## con1       0.80      0.81    0.78      0.46 4.2    0.017 0.0037  0.46
## con3       0.79      0.80    0.77      0.45 4.0    0.018 0.0018  0.45
## con6       0.79      0.80    0.77      0.45 4.0    0.018 0.0032  0.45
## con8       0.79      0.80    0.77      0.44 4.0    0.018 0.0039  0.43
## con9       0.78      0.79    0.76      0.44 3.9    0.018 0.0022  0.45
## con10      0.82      0.82    0.79      0.48 4.6    0.015 0.0016  0.48
## 
##  Item statistics 
##         n raw.r std.r r.cor r.drop mean   sd
## con1  343  0.70  0.72  0.64   0.58  4.0 0.81
## con3  343  0.74  0.75  0.68   0.61  3.9 0.92
## con6  343  0.74  0.75  0.68   0.61  3.8 0.91
## con8  343  0.74  0.76  0.69   0.63  4.0 0.83
## con9  343  0.77  0.77  0.72   0.65  3.7 0.96
## con10 343  0.71  0.68  0.57   0.52  3.1 1.17
## 
## Non missing response frequency for each item
##          1    2    3    4    5 miss
## con1  0.00 0.02 0.26 0.42 0.30    0
## con3  0.01 0.06 0.26 0.37 0.31    0
## con6  0.00 0.07 0.31 0.36 0.27    0
## con8  0.00 0.04 0.22 0.44 0.30    0
## con9  0.01 0.09 0.29 0.37 0.24    0
## con10 0.10 0.22 0.29 0.26 0.12    0

Check raw_alpha for the overall reliability. In the Item Statistics table, look at the Drop column — if any value exceeds the overall alpha, dropping that item would improve reliability.

Exercise 3

pe_cs <- parameterEstimates(fit_con_sup, standardized = TRUE)
loadings <- pe_cs[pe_cs$op == "=~", c("lhs", "rhs", "est", "std.all")]
loadings
##    lhs   rhs   est std.all
## 1  con  con1 1.000   0.646
## 2  con  con3 1.227   0.693
## 3  con  con6 1.208   0.692
## 4  con  con8 1.110   0.697
## 5  con  con9 1.333   0.725
## 6  con con10 1.289   0.576
## 7  sup  sup2 1.000   0.670
## 8  sup  sup4 1.290   0.726
## 9  sup  sup5 1.023   0.701
## 10 sup  sup7 0.907   0.520

For each factor, the row with the largest std.all value is the strongest indicator.

Exercise 4

modindices(fit_con_sup, minimum.value = 3.84, sort. = TRUE)
##     lhs op  rhs    mi    epc sepc.lv sepc.all sepc.nox
## 58 con8 ~~ con9 9.258 -0.089  -0.089   -0.227   -0.227
## 45 con3 ~~ con9 7.828  0.091   0.091    0.207    0.207

If any rows appear, examine the lhs, op, and rhs columns to see which parameter would be freed. Theoretical justification depends on the substance of the items — for example, freeing a residual covariance between two items with similar wording could be justified by a common-method interpretation.

Summary & next tutorial

You have now fitted a complete CFA, read every section of the output, extracted fit indices and parameter estimates, computed Cronbach’s alpha, examined modification indices, drawn a path diagram, and written up the results in APA format. This is the workflow you will repeat for every CFA in your dissertation. In the next tutorial we extend the CFA toolkit: we build parcels from individual items and fit a second-order factor model, both of which appear frequently in published organisational psychology research.

Further reading

  • Rosseel (2012) — the lavaan reference paper, covering the full range of models the package can fit.
  • Beaujean (2014) — a lavaan-based walkthrough of CFA, including parceling and second-order factors.
  • Kelloway (2015) — the textbook on which the original Mplus teaching materials were based; useful for cross-referencing Mplus output with lavaan output.
  • Hu & Bentler (1999) — the source of the fit-index cutoff criteria used in the write-up.

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