`compare_performance()`: Is there a way to get robust fit indices of lavaan models fitted using robust maximum likelihood?

[HanZhang-psych]

lavaan models fitted using robust maximum likelihood give the following fit indices:

compare_performance() appears to use the fit indices under maximum likelihood (e.g., for CFI, .874).

Is there a way to get robust fit indices of lavaan models fitted using robust maximum likelihood (i.e., .875 in the picture)?

Thanks!

[bwiernik]

Can you give a reproducible example for a model and getting robust indices with lavaan?

[HanZhang-psych]

Sure!

library(lavaan)
library(performance)

model <- "
 visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
  speed =~ x7 + x8 + x9"

fit <- cfa(model, data=HolzingerSwineford1939, estimator='MLM')

Model summary:

> summary(fit, fit.measures=T)
lavaan 0.6-8 ended normally after 35 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        21

  Number of observations                           301

Model Test User Model:
                                              Standard      Robust
  Test Statistic                                85.306      80.872
  Degrees of freedom                                24          24
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.055
       Satorra-Bentler correction                                 

Model Test Baseline Model:

  Test statistic                               918.852     789.298
  Degrees of freedom                                36          36
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.164

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.931       0.925
  Tucker-Lewis Index (TLI)                       0.896       0.887

  Robust Comparative Fit Index (CFI)                         0.932
  Robust Tucker-Lewis Index (TLI)                            0.897

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -3737.745   -3737.745
  Loglikelihood unrestricted model (H1)      -3695.092   -3695.092

  Akaike (AIC)                                7517.490    7517.490
  Bayesian (BIC)                              7595.339    7595.339
  Sample-size adjusted Bayesian (BIC)         7528.739    7528.739

Root Mean Square Error of Approximation:

  RMSEA                                          0.092       0.089
  90 Percent confidence interval - lower         0.071       0.068
  90 Percent confidence interval - upper         0.114       0.110
  P-value RMSEA <= 0.05                          0.001       0.001

  Robust RMSEA                                               0.091
  90 Percent confidence interval - lower                     0.070
  90 Percent confidence interval - upper                     0.113

Standardized Root Mean Square Residual:

  SRMR                                           0.065       0.065

Parameter Estimates:

  Standard errors                           Robust.sem
  Information                                 Expected
  Information saturated (h1) model          Structured

Alternatively:

> fitMeasures(fit, fit.measures = c('chisq','cfi','tli','rmsea')) # ML (standard) estimates
 chisq    cfi    tli  rmsea 
85.306  0.931  0.896  0.092 

> fitMeasures(fit, fit.measures = c('chisq.scaled','cfi.robust','tli.robust','rmsea.robust')) # MLM (robust) estimates
chisq.scaled   cfi.robust   tli.robust rmsea.robust 
      80.872        0.932        0.897        0.091 

compare_performance() seems to only give standard estimates:

> compare_performance(fit, metrics='all')
# Comparison of Model Performance Indices

Name |  Model | Chi2(24) | p (Chi2) | Baseline(36) | p (Baseline) |   GFI |  AGFI |   NFI |  NNFI |   CFI | RMSEA |    RMSEA  CI | p (RMSEA) |   RMR |  SRMR |   RFI |  PNFI |   IFI |   RNI | Loglikelihood |      AIC | AIC weights |      BIC | BIC weights | BIC_adjusted
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
fit  | lavaan |   85.306 |   < .001 |      918.852 |       < .001 | 0.943 | 0.894 | 0.907 | 0.896 | 0.931 | 0.092 | [0.07, 0.11] |    < .001 | 0.082 | 0.065 | 0.861 | 0.605 | 0.931 | 0.931 |     -3737.745 | 7517.490 |        1.00 | 7595.339 |        1.00 |     7528.739

However, reporting the robust fit indices is desirable when data is not multivariate normal.

Also, FWIW, when using MLM, lavaan also gives you "scaled" alternative fit indices:

> fitMeasures(fit, fit.measures = c('cfi.scaled','tli.scaled','rmsea.scaled'))
  cfi.scaled   tli.scaled rmsea.scaled 
       0.925        0.887        0.089 

But it seems like the "robust" versions are more desirable (see: https://lavaan.ugent.be/history/dot5.html, Version 0.5-21)

[DominiqueMakowski]

that seems like a good idea to add. @HanZhang-psych would you try making a PR? The code is pretty straightforward:

https://github.com/easystats/performance/blob/b181ee03eaa78075825d79dda725ce8eaf928ae3/R/model_performance.lavaan.R#L114-L117

We could also make a metrics="robust" option to get only these

[bwiernik]

What is you specify method="MLR"? I think that's the main way you ask for these in lavaan? At estimation time, not at summary()

[HanZhang-psych]

@bwiernik yes, the estimator is specified when fitting the model. e.g., fit <- cfa(model, data=HolzingerSwineford1939, estimator='MLM'). lavaan provides several 'robust' estimators, including MLM and MLR.