Robust CFI Discrepancy

[Garett Howardson]

Good morning, 

I'm running a full SEM correcting for non-normality with Satorra-Bentler. In trying to confirm some of the robust statistics, I keep running into an issue in that the robust CFI I'm calculating using information in the lavaan "test" slot differs from the value reported by the print.fit.measures() function as seen here: https://github.com/yrosseel/lavaan/blob/master/R/lav_fit_measures.R. 

From the lav_fit_measures() function shown at the link above, the robust CFI is calculated (see code line 373) using using Equation (15) here. 

The subroutine for calculating the robust CFI is given as: 

if("cfi.robust" %in% fit.measures) {
	if(TEST[[2]]$test %in% c("satorra.bentler", "yuan.bentler")) {
	# see Brosseau-Liard & Savalei MBR 2014, equation 15
	# what to do if X2 = 0 and df = 0? in this case,
	# the scaling factor (ch) will be NA, and we get NA
	# (instead of 1)
	if(X2 < .Machine$double.eps && df == 0) {
	ch <- 0
	} else {
	ch <- TEST[[2]]$scaling.factor
	}
	cb <- fit.indep@test[[2]]$scaling.factor
	t1 <- max( c(X2 - (ch*df), 0) )
	t2 <- max( c(X2 - (ch*df), X2.null - (cb*df.null), 0) )
	if(is.na(t1) || is.na(t2)) {
	indices["cfi.robust"] <- NA
	} else if(t1 == 0 && t2 == 0) {
	indices["cfi.robust"] <- 1
	} else {
	indices["cfi.robust"] <- 1 - t1/t2
	}
	} else {
	indices["cfi.robust"] <- NA
	}
	}

I'm looking specifically at t1 and t2 values. The equation given in Brosseau-Liard & Savalei MBR 2014, equation 15 is

I'm reading this as T_ML,H is the Satorra-Bentler scaled chi-square or T statistic for the hypothesized model and the T_ML,B is the same quantity for the baseline / null model. the c quantities are the respective correction factors for the theoretical and baseline / null models, respectively and df is likewise degrees of freedom. 

If my fitted hypothesis model object were named epm.sem.fit.d.mlm, I'm drawing the said quantities from the epm.sem.fit.d.mlm@test[[2]] slot of the fitted object with the T value coming from the epm.sem.fit.d.mlm@test[[2]]$stat, degrees of freedom from epm.sem.fit.d.mlm@test[[2]]$df, and Satorra-Bentler scaling factor from epm.sem.fit.d.mlm@test[[2]]$scaling.factor. In other words,

chisq.fit = epm.sem.fit.d.mlm@test[[2]]$stat
df.fit = epm.sem.fit.d.mlm@test[[2]]$df
sf.fit = epm.sem.fit.d.mlm@test[[2]]$scaling.factor

If the fitted null (independence model) object were named epm.ov.null.fit, I'm drawing the relevant quantities from epm.ov.null.fit@test[[2]]$stat, epm.ov.null.fit@test[[2]]$df, and epm.ov.null.fit@test[[2]]$scaling.factor. In other words, 

chisq.null = epm.ov.null.fit@test[[2]]$stat
df.null = epm.ov.null.fit@test[[2]]$df
sf.null = epm.ov.null.fit@test[[2]]$scaling.factor

If the numerator and denominator of the equation above were named cfi.num and cfi.denom, respectively, I would calculate these quantities as

cfi.num = (chisq.fit - df.fit * sf.fit)
cfi.denom = (chisq.null - df.null * sf.null)

The robust CFI (cfi.rob) value would then be

cfi.rob = 1 - (cfi.num / cfi.denom)

It looks like the function calculating this value looks for the maximum between the cfi.num and cfi.denom and chooses that as the actual denominator to avoid producing values beyond 1, which makes perfect sense. The actual values associated with each of the quantities above when drawing from the fitted object are:

chisq.fit
[1] 1826.315

> df.fit
[1] 1105

> sf.fit
[1] 1.258276

> chisq.null
[1] 5791.451

> df.null
[1] 1176

> sf.null
[1] 1.390944

> cfi.num
[1] 435.9202

> cfi.denom
[1] 4155.701

Plugging this information into the above then results in a robust cfi value of 

cfi.rob
[1] 0.8951031

or approximately .90 when rounded. When calling the summary() function and passing in the fitted lavaan object, which I believe calls the print.fit.measures() function noted above, the snippet of results relevant to the issue returned is 

summary(epm.sem.fit.d.mlm, fit.measures = TRUE)
lavaan (0.5-23.1097) converged normally after 101 iterations


                                                  Used       Total
  Number of observations                           231         254


  Estimator                                         ML      Robust
  Minimum Function Test Statistic             2298.008    1826.315
  Degrees of freedom                              1105        1105
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.258
    for the Satorra-Bentler correction


Model test baseline model:


  Minimum Function Test Statistic             8055.582    5791.443
  Degrees of freedom                              1176        1176
  P-value                                        0.000       0.000


User model versus baseline model:


  Comparative Fit Index (CFI)                    0.827       0.844
  Tucker-Lewis Index (TLI)                       0.815       0.834


  Robust Comparative Fit Index (CFI)                         0.859
  Robust Tucker-Lewis Index (TLI)                            0.850

The Robust CFI values reported by the summary call passing in the lavaan fit object prints a robust value of .859, which seems to have transposed the second and third decimal. I should note that the print null chi square is slightly different than the value extracted from running the null model myself, but the differences are minor. Using the printed result instead of the extracted result still produces a robust CFI value of .895. I've also tried calling fitMeasures() directly instead of the summary() function and the results returned match the results presented by the summary() function: 

fmsr = fitMeasures(epm.sem.fit.d.mlm)
as.data.frame(cbind(fit.meas = names(fmsr), val = fmsr), row.names = FALSE)
                        fit.meas                val
1                           npar                169
2                           fmin   4.97404358398868
3                          chisq   2298.00813580277
4                             df               1105
5                         pvalue                  0
6                   chisq.scaled    1826.3150123355
7                      df.scaled               1105
8                  pvalue.scaled                  0
9           chisq.scaling.factor   1.25827588355859
10                baseline.chisq   8055.58155142839
11                   baseline.df               1176
12               baseline.pvalue                  0
13         baseline.chisq.scaled   5791.44276380469
14            baseline.df.scaled               1176
15        baseline.pvalue.scaled                  0
16 baseline.chisq.scaling.factor   1.39094555190532
17                           cfi  0.826587107531988
18                           tli  0.815444740685627
19                          nnfi  0.815444740685627
20                           rfi  0.696401440220714
21                           nfi  0.714730945105348
22                          pnfi  0.671579672059022
23                           ifi  0.828358515474493
24                           rni  0.826587107531988
25                    cfi.scaled  0.843717049642081
26                    tli.scaled  0.833675339709582
27                    cfi.robust  0.858623461445068
28                    tli.robust  0.849539539058281
29                   nnfi.scaled  0.833675339709582
30                   nnfi.robust  0.849539539058281
31                    rfi.scaled  0.664390717641696
32                    nfi.scaled  0.684652842682036
33                    ifi.scaled  0.684652842682036
34                    rni.scaled  0.895151324692744
35                    rni.robust  0.858623461445068
36                          logl  -12516.4824760447
37             unrestricted.logl  -11367.4784081433
38                           aic   25370.9649520894
39                           bic   25952.7335451676
40                        ntotal                231
41                          bic2   25417.0993534133
42                         rmsea 0.0683651189435794
43                rmsea.ci.lower  0.064424371943673
44                rmsea.ci.upper 0.0723016042892283
45                  rmsea.pvalue                  0
46                  rmsea.scaled 0.0531588229969579
47         rmsea.ci.lower.scaled 0.0492753119101313
48         rmsea.ci.upper.scaled 0.0569936248456142
49           rmsea.pvalue.scaled 0.0895723646043242
50                  rmsea.robust 0.0596297918004684
51         rmsea.ci.lower.robust 0.0547387450182823
52         rmsea.ci.upper.robust 0.0644519622543453
53           rmsea.pvalue.robust               <NA>
54                           rmr 0.0696247719099917
55                    rmr_nomean 0.0696247719099917
56                          srmr 0.0723072013115454
57                  srmr_bentler 0.0723072013115454
58           srmr_bentler_nomean 0.0737391660921703
59                   srmr_bollen 0.0722919658035523
60            srmr_bollen_nomean 0.0736989433759451
61                    srmr_mplus 0.0722967647188455
62             srmr_mplus_nomean 0.0737038389676122
63                         cn_05   119.962150339531
64                         cn_01   123.364744748049
65                           gfi  0.953218803530716
66                          agfi  0.946064032306002
67                          pgfi  0.826771411225621
68                           mfi 0.0756022909855766
69                          ecvi               <NA>

Any thoughts on why the values reported by fitMeasures() and summary() differ from those calculated directly from the null and fitted objects' slots? I have a sneaking suspicion that I'm missing something incredibly simple, but I just want to confirm before moving forward with the paper. Please let me know if you need any additional information. 

Thank you!

Garett

[Terrence Jorgensen]

I'm reading this as T_ML,H is the Satorra-Bentler scaled chi-square or T statistic for the hypothesized model and the T_ML,B is the same quantity for the baseline / null model. 

I think this is your problem.  The ML subscript means it is the naïve chi-squared.  Notice in eq. 11 that the subscript SB is used to denote the scaled test statistic.

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

UvA web page: http://www.uva.nl/profile/t.d.jorgensen

[Yves Rosseel]

I agree with Terrence. T_ML,H is NOT the Satorra-Bentler scaled test
statistic. It is the standard test statistic, which you need to extract
from the @test[[1]] slot. As far as I can tell, the implementation in
lavaan is correct.

Nevertheless, I appreciate it enormously that you took the effort to
manually check these computations. This is how we ensure code quality in
the open-source world.

Yves.

On 11/08/2017 11:22 AM, Terrence Jorgensen wrote:
>
> <https://lh3.googleusercontent.com/-FdN0bqOMf94/WgGlK3LA4SI/AAAAAAAADAM/tl8RfWSeJcsCBSDFd2PQmf9R271bey5cACLcBGAs/s1600/Screen%2BShot%2B2017-11-07%2Bat%2B7.20.26%2BAM.png>

> --
> You received this message because you are subscribed to the Google
> Groups "lavaan" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to lavaan+un...@googlegroups.com
> <mailto:lavaan+un...@googlegroups.com>.
> To post to this group, send email to lav...@googlegroups.com
> <mailto:lav...@googlegroups.com>.
> Visit this group at https://groups.google.com/group/lavaan.
> For more options, visit https://groups.google.com/d/optout.

[Garett Howardson]

Yep. I definitely missed something “incredibly simple.” Thank you both for your help - it is MUCH appreciated.