R/lav_fit_cfi.R

In lavaan: Latent Variable Analysis

# functions related to CFI and other 'incremental' fit indices

# lower-level functions:
# - lav_fit_cfi
# - lav_fit_rni (same as CFI, but without the max(0,))
# - lav_fit_tli/lav_fit_nnfi
# - lav_fit_rfi
# - lav_fit_nfi
# - lav_fit_pnfi
# - lav_fit_ifi

# higher-level functions:
# - lav_fit_cfi_lavobject

# Y.R. 20 July 2022

# CFI - comparative fit index (Bentler, 1990)
# robust version: Brosseau-Liard & Savalei MBR 2014, equation 15

# robust version MLMV (scaled.shifted)
# Savalei, V. (2018). On the computation of the RMSEA and CFI from the
# mean-and-variance corrected test statistic with nonnormal data in SEM.
# Multivariate behavioral research, 53(3), 419-429. eq 9

# note: robust MLM == robust MLMV

# when missing = "fiml":
# Zhang, X., & Savalei, V. (2023). New computations for RMSEA and CFI following
# FIML and TS estimation with missing data. Psychological Methods, 28(2),
# 263-283. https://doi.org/10.1037/met0000445


lav_fit_cfi <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL,
                        c.hat = 1, c.hat.null = 1) {

    if(anyNA(c(X2, df, X2.null, df.null, c.hat, c.hat.null))) {
        return(as.numeric(NA))
    }

    # robust?
    if(df > 0 && !missing(c.hat) && !missing(c.hat.null) &&
       c.hat != 1 && c.hat.null != 1) {
        t1 <- max( c(X2 - (c.hat * df), 0) )
        t2 <- max( c(X2 - (c.hat * df), X2.null - (c.hat.null * df.null), 0) )
    } else {
        t1 <- max( c(X2 - df, 0) )
        t2 <- max( c(X2 - df, X2.null - df.null, 0) )
    }

    if(isTRUE(all.equal(t1, 0)) && isTRUE(all.equal(t2, 0))) {
        CFI <- 1
    } else {
        CFI <- 1 - t1/t2
    }

    CFI
}

# RNI - relative noncentrality index (McDonald & Marsh, 1990)
# same as CFI, but without the max(0,)
lav_fit_rni <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL,
                        c.hat = 1, c.hat.null = 1) {

    if(anyNA(c(X2, df, X2.null, df.null, c.hat, c.hat.null))) {
        return(as.numeric(NA))
    }

    # robust?
    if(df > 0 && !missing(c.hat) && !missing(c.hat.null) &&
       c.hat != 1 && c.hat.null != 1) {
        t1 <- X2 - (c.hat * df)
        t2 <- X2.null - (c.hat.null * df.null)
    } else {
        t1 <- X2 - df
        t2 <- X2.null - df.null
    }

    if(isTRUE(all.equal(t2, 0))) {
        RNI <- as.numeric(NA)
    } else if(!is.finite(t1) || !is.finite(t2)) {
        RNI <- as.numeric(NA)
    } else {
        RNI <- 1 - t1/t2
    }

    RNI
}

# TLI - Tucker-Lewis index (Tucker & Lewis, 1973)
# same as
# NNFI - nonnormed fit index (NNFI, Bentler & Bonett, 1980)
# note: formula in lavaan <= 0.5-20:
# t1 <- X2.null/df.null - X2/df
# t2 <- X2.null/df.null - 1
# if(t1 < 0 && t2 < 0) {
#    TLI <- 1
#} else {
#    TLI <- t1/t2
#}
# note: TLI original formula was in terms of fx/df, not X2/df
# then, t1 <- fx_0/df.null - fx/df
#       t2 <- fx_0/df.null - 1/N (or N-1 for wishart)

# note: in lavaan 0.5-21, we use the alternative formula:
# TLI <- 1 - ((X2 - df)/(X2.null - df.null) * df.null/df)
# - this one has the advantage that a 'robust' version
#   can be derived; this seems non-trivial for the original one
# - unlike cfi, we do not use 'max(0, )' for t1 and t2
#   therefore, t1 can go negative, and TLI can be > 1
lav_fit_tli <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL,
                        c.hat = 1, c.hat.null = 1) {

    if(anyNA(c(X2, df, X2.null, df.null, c.hat, c.hat.null))) {
        return(as.numeric(NA))
    }

    # robust?
    if(df > 0 && !missing(c.hat) && !missing(c.hat.null) &&
       c.hat != 1 && c.hat.null != 1) {
        t1 <- (X2      - c.hat      * df     ) * df.null
        t2 <- (X2.null - c.hat.null * df.null) * df
    } else {
        t1 <- (X2      - df)      * df.null
        t2 <- (X2.null - df.null) * df
    }

    if(df > 0 && abs(t2) > 0) {
        TLI <- 1 - t1/t2
    } else if(!is.finite(t1) || !is.finite(t2)) {
        TLI <- as.numeric(NA)
    } else {
       TLI <- 1
    }

    TLI
}

# alias for nnfi
lav_fit_nnfi <- lav_fit_tli

# RFI - relative fit index (Bollen, 1986; Joreskog & Sorbom 1993)
lav_fit_rfi <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL) {

    if(anyNA(c(X2, df, X2.null, df.null))) {
        return(as.numeric(NA))
    }

    if(df > df.null) {
        RLI <- as.numeric(NA)
    } else if(df > 0 && df.null > 0) {
        t1 <- X2.null/df.null - X2/df
        t2 <- X2.null/df.null
        if(!is.finite(t1) || !is.finite(t2)) {
            RLI <- as.numeric(NA)
        } else if(t1 < 0 || t2 < 0) {
            RLI <- 1
        } else {
            RLI <- t1/t2
        }
    } else {
       RLI <- 1
    }

    RLI
}

# NFI - normed fit index (Bentler & Bonett, 1980)
lav_fit_nfi <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL) {

    if(anyNA(c(X2, df, X2.null, df.null))) {
        return(as.numeric(NA))
    }

    if(df > df.null || isTRUE(all.equal(X2.null,0))) {
        NFI <- as.numeric(NA)
    } else if(df > 0) {
        t1 <- X2.null - X2
        t2 <- X2.null
        NFI <- t1/t2
    } else {
        NFI <- 1
    }

    NFI
}

# PNFI - Parsimony normed fit index (James, Mulaik & Brett, 1982)
lav_fit_pnfi <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL) {

    if(anyNA(c(X2, df, X2.null, df.null))) {
        return(as.numeric(NA))
    }

    if(df.null > 0 && X2.null > 0) {
        t1 <- X2.null - X2
        t2 <- X2.null
        PNFI <- (df/df.null) * (t1/t2)
    } else {
        PNFI <- as.numeric(NA)
    }

    PNFI
}

# IFI - incremental fit index (Bollen, 1989; Joreskog & Sorbom, 1993)
lav_fit_ifi <- function(X2 = NULL, df = NULL, X2.null = NULL, df.null = NULL) {

    if(anyNA(c(X2, df, X2.null, df.null))) {
        return(as.numeric(NA))
    }

    t1 <- X2.null - X2
    t2 <- X2.null - df
    if(!is.finite(t1) || !is.finite(t2)) {
        IFI <- as.numeric(NA)
    } else if(t2 < 0) {
        IFI <- 1
    } else if(isTRUE(all.equal(t2,0))) {
        IFI <- as.numeric(NA)
    } else {
        IFI <- t1/t2
    }

    IFI
}

# higher-level function
lav_fit_cfi_lavobject <- function(lavobject = NULL, fit.measures = "cfi",
                                  baseline.model = NULL,
                                  standard.test = "standard",
                                  scaled.test   = "none",
                                  robust        = TRUE,
                                  cat.check.pd  = TRUE) {

    # check lavobject
    stopifnot(inherits(lavobject, "lavaan"))

    # check for categorical
    categorical.flag <- lavobject@Model@categorical

    # tests
    TEST <- lavobject@test
    test.names <- sapply(lavobject@test, "[[", "test")
    if(test.names[1] == "none" || standard.test == "none") {
        return(list())
    }
    test.idx <- which(test.names == standard.test)[1]
    if(length(test.idx) == 0L) {
        return(list())
    }

    scaled.flag <- FALSE
    if(!scaled.test %in% c("none", "standard", "default")) {
        scaled.idx <- which(test.names == scaled.test)
        if(length(scaled.idx) > 0L) {
            scaled.idx <- scaled.idx[1] # only the first one
            scaled.flag <- TRUE
        }
    }

    # robust?
    robust.flag <- FALSE
    if(robust && scaled.flag &&
       scaled.test %in% c("satorra.bentler", "yuan.bentler.mplus",
                          "yuan.bentler", "scaled.shifted")) {
        robust.flag <- TRUE
    }

    # FIML?
    fiml.flag <- FALSE
    if(robust && lavobject@Options$missing %in% c("ml", "ml.x")) {
        fiml.flag <- robust.flag <- TRUE
        # check if we can compute corrected values
        if(scaled.flag) {
            version <- "V3"
        } else {
            version <- "V6"
        }
        fiml <- try(lav_fit_fiml_corrected(lavobject, version = version),
                    silent = TRUE)
        if(inherits(fiml, "try-error")) {
            warning("lavaan WARNING: computation of robust CFI failed.")
            fiml <- list(XX3 = as.numeric(NA), df3 = as.numeric(NA),
                         c.hat3= as.numeric(NA), XX3.scaled = as.numeric(NA),
                         XX3.null = as.numeric(NA), df3.null = as.numeric(NA),
                         c.hat3.null = as.numeric(NA))
        } else if(anyNA(c(fiml$XX3, fiml$df3, fiml$c.hat3, fiml$XX3.scaled,
                    fiml$XX3.null, fiml$df3.null, fiml$c.hat3.null))) {
            warning("lavaan WARNING: computation of robust CFI resulted in NA values.")
        }
    }

    # supported fit measures in this function
    # baseline model
    fit.baseline <- c("baseline.chisq", "baseline.df", "baseline.pvalue")
    if(scaled.flag) {
        fit.baseline <- c(fit.baseline, "baseline.chisq.scaled",
                          "baseline.df.scaled", "baseline.pvalue.scaled",
                          "baseline.chisq.scaling.factor")
    }

    fit.cfi.tli <- c("cfi", "tli")
    if(scaled.flag) {
        fit.cfi.tli <- c(fit.cfi.tli, "cfi.scaled", "tli.scaled")
    }
    if(robust.flag) {
        fit.cfi.tli <- c(fit.cfi.tli, "cfi.robust", "tli.robust")
    }

    # other incremental fit indices
    fit.cfi.other <- c("nnfi", "rfi", "nfi", "pnfi", "ifi", "rni")
    if(scaled.flag) {
        fit.cfi.other <- c(fit.cfi.other, "nnfi.scaled", "rfi.scaled",
                       "nfi.scaled", "pnfi.scaled", "ifi.scaled", "rni.scaled")
    }
    if(robust.flag) {
        fit.cfi.other <- c(fit.cfi.other, "nnfi.robust", "rni.robust")
    }

    # which one do we need?
    if(missing(fit.measures)) {
        # default set
        fit.measures <- c(fit.baseline, fit.cfi.tli)
    } else {
        # remove any not-CFI related index from fit.measures
        rm.idx <- which(!fit.measures %in%
                        c(fit.baseline, fit.cfi.tli, fit.cfi.other))
        if(length(rm.idx) > 0L) {
            fit.measures <- fit.measures[-rm.idx]
        }
        if(length(fit.measures) == 0L) {
            return(list())
        }
    }


    # basic test statistics
    X2 <- TEST[[test.idx]]$stat
    df <- TEST[[test.idx]]$df
    G <- lavobject@Data@ngroups  # number of groups
    N <- lav_utils_get_ntotal(lavobject = lavobject) # N vs N-1

    # scaled X2
    if(scaled.flag) {
        X2.scaled <- TEST[[scaled.idx]]$stat
        df.scaled <- TEST[[scaled.idx]]$df
    }
    if(robust.flag) {
        XX3 <- X2
        if(categorical.flag) {
            out <- try(lav_fit_catml_dwls(lavobject, check.pd = cat.check.pd),
                       silent = TRUE)
            if(inherits(out, "try-error")) {
                XX3 <- df3 <- c.hat <- c.hat3 <- XX3.scaled <- as.numeric(NA)
            } else {
                XX3 <- out$XX3
                df3 <- out$df3
                c.hat3 <- c.hat <- out$c.hat3
                XX3.scaled <- out$XX3.scaled
            }
        } else if(fiml.flag) {
            XX3 <- fiml$XX3
            df3 <- fiml$df3
            c.hat3 <- c.hat <- fiml$c.hat3
            XX3.scaled <- fiml$XX3.scaled
        } else if(scaled.test == "scaled.shifted") {
            # compute c.hat from a and b
            a <- TEST[[scaled.idx]]$scaling.factor
            b <- TEST[[scaled.idx]]$shift.parameter
            c.hat <- a * (df - b) / df
        } else {
            c.hat <- TEST[[scaled.idx]]$scaling.factor
        }
    }

    # output container
    indices <- list()

    # only do what is needed (per groups)
    cfi.baseline.flag <- cfi.tli.flag <- cfi.other.flag <- FALSE
    if(any(fit.baseline %in% fit.measures)) {
        cfi.baseline.flag <- TRUE
    }
    if(any(fit.cfi.tli %in% fit.measures)) {
        cfi.tli.flag <- TRUE
    }
    if(any(fit.cfi.other %in% fit.measures)) {
        cfi.other.flag <- TRUE
    }

    # 1. BASELINE model
    baseline.test <- NULL

    # we use the following priority:
    # 1. user-provided baseline model
    # 2. baseline model in @external slot
    # 3. baseline model in @baseline slot
    # 4. nothing -> compute independence model

    # 1. user-provided baseline model
    if( !is.null(baseline.model) ) {
        baseline.test <-
            lav_fit_measures_check_baseline(fit.indep = baseline.model,
                                            object    = lavobject)
    # 2. baseline model in @external slot
    } else if( !is.null(lavobject@external$baseline.model) ) {
        fit.indep <- lavobject@external$baseline.model
        baseline.test <-
            lav_fit_measures_check_baseline(fit.indep = fit.indep,
                                            object    = lavobject)
    # 3. internal @baseline slot
    } else if( .hasSlot(lavobject, "baseline") &&
               length(lavobject@baseline) > 0L &&
               !is.null(lavobject@baseline$test) ) {
        baseline.test <- lavobject@baseline$test
    # 4. (re)compute independence model
    } else {
        fit.indep <- try(lav_object_independence(lavobject), silent = TRUE)
        baseline.test <-
            lav_fit_measures_check_baseline(fit.indep = fit.indep,
                                            object    = lavobject)
    }

    # baseline.test.idx
    baseline.test.idx <- which(names(baseline.test) == standard.test)[1]
    if(scaled.flag) {
        baseline.scaled.idx <- which(names(baseline.test) == scaled.test)[1]
    }

    if(!is.null(baseline.test)) {
        X2.null <- baseline.test[[baseline.test.idx]]$stat
        df.null <- baseline.test[[baseline.test.idx]]$df
        if(scaled.flag) {
            X2.null.scaled <- baseline.test[[baseline.scaled.idx]]$stat
            df.null.scaled <- baseline.test[[baseline.scaled.idx]]$df
        }
        if(robust.flag) {
            XX3.null <- X2.null
            if(categorical.flag) {
                if(inherits(out, "try-error")) {
                    XX3.null <- c.hat.null <- as.numeric(NA)
                } else {
                    XX3.null   <- out$XX3.null
                    c.hat.null <- out$c.hat3.null
                }
            } else if(fiml.flag) {
                XX3.null   <- fiml$XX3.null
                c.hat.null <- fiml$c.hat3.null
            } else if(scaled.test == "scaled.shifted") {
                # compute c.hat from a and b
                a.null <-
                    baseline.test[[baseline.scaled.idx]]$scaling.factor
                b.null <-
                    baseline.test[[baseline.scaled.idx]]$shift.parameter
                c.hat.null <- a.null * (df.null - b.null) / df.null
            } else {
                c.hat.null <-
                    baseline.test[[baseline.scaled.idx]]$scaling.factor
            }
        }
    } else {
        X2.null <- df.null <- as.numeric(NA)
        X2.null.scaled <- df.null.scaled <- as.numeric(NA)
        c.hat.null <- as.numeric(NA)
    }

    # check for NAs of nonfinite numbers
    if(!is.finite(X2) || !is.finite(df) ||
       !is.finite(X2.null) || !is.finite(df.null)) {
        indices[fit.measures] <- as.numeric(NA)
        return(indices)
    }

    # fill in baseline indices
    if(cfi.baseline.flag) {
        indices["baseline.chisq"]  <- X2.null
        indices["baseline.df"]     <- df.null
        indices["baseline.pvalue"] <- baseline.test[[baseline.test.idx]]$pvalue
        if(scaled.flag) {
            indices["baseline.chisq.scaled"]  <- X2.null.scaled
            indices["baseline.df.scaled"]     <- df.null.scaled
            indices["baseline.pvalue.scaled"] <-
                baseline.test[[baseline.scaled.idx]]$pvalue
            indices["baseline.chisq.scaling.factor"] <-
                baseline.test[[baseline.scaled.idx]]$scaling.factor
        }
    }

    # 2. CFI and TLI
    if(cfi.tli.flag) {
        indices["cfi"] <- lav_fit_cfi(X2 = X2, df = df,
                                      X2.null = X2.null, df.null = df.null)
        indices["tli"] <- lav_fit_tli(X2 = X2, df = df,
                                      X2.null = X2.null, df.null = df.null)
        if(scaled.flag) {
            indices["cfi.scaled"] <-
                lav_fit_cfi(X2 = X2.scaled, df = df.scaled,
                            X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["tli.scaled"] <-
                lav_fit_tli(X2 = X2.scaled, df = df.scaled,
                            X2.null = X2.null.scaled, df.null = df.null.scaled)
        }
        if(robust.flag) {
            indices["cfi.robust"] <-
                lav_fit_cfi(X2 = XX3, df = df,
                            X2.null = XX3.null, df.null = df.null,
                            c.hat = c.hat, c.hat.null = c.hat.null)
            indices["tli.robust"] <-
                lav_fit_tli(X2 = XX3, df = df,
                            X2.null = XX3.null, df.null = df.null,
                            c.hat = c.hat, c.hat.null = c.hat.null)
        }
    }

    # 3. other
    # c("nnfi", "rfi", "nfi", "pnfi", "ifi", "rni")
    if(cfi.other.flag) {
        indices["nnfi"] <-
            lav_fit_nnfi(X2 = X2, df = df, X2.null = X2.null, df.null = df.null)
        indices["rfi"] <-
            lav_fit_rfi( X2 = X2, df = df, X2.null = X2.null, df.null = df.null)
        indices["nfi"] <-
            lav_fit_nfi( X2 = X2, df = df, X2.null = X2.null, df.null = df.null)
        indices["pnfi"] <-
            lav_fit_pnfi(X2 = X2, df = df, X2.null = X2.null, df.null = df.null)
        indices["ifi"] <-
            lav_fit_ifi( X2 = X2, df = df, X2.null = X2.null, df.null = df.null)
        indices["rni"] <-
            lav_fit_rni( X2 = X2, df = df, X2.null = X2.null, df.null = df.null)

        if(scaled.flag) {
            indices["nnfi.scaled"] <-
                lav_fit_nnfi(X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["rfi.scaled"] <-
                lav_fit_rfi( X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["nfi.scaled"] <-
                lav_fit_nfi( X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["pnfi.scaled"] <-
                lav_fit_pnfi(X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["ifi.scaled"] <-
                lav_fit_ifi( X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
            indices["rni.scaled"] <-
                lav_fit_rni( X2 = X2.scaled, df = df.scaled,
                             X2.null = X2.null.scaled, df.null = df.null.scaled)
        }
        if(robust.flag) {
            indices["nnfi.robust"] <-
                lav_fit_nnfi(X2 = XX3, df = df,
                             X2.null = XX3.null, df.null = df.null,
                             c.hat = c.hat, c.hat.null = c.hat.null)
            indices["rni.robust"] <-
                lav_fit_rni(X2 = XX3, df = df,
                            X2.null = XX3.null, df.null = df.null,
                            c.hat = c.hat, c.hat.null = c.hat.null)
        }
    }

    # return only those that were requested
    indices[fit.measures]
}


# new in 0.6-5
# internal function to check the (external) baseline model, and
# return baseline 'test' list if everything checks out (and NULL otherwise)
lav_fit_measures_check_baseline <- function(fit.indep = NULL, object = NULL) {

    TEST <- NULL

    # check if everything is in order
    if( inherits(fit.indep, "try-error") ) {
        warning("lavaan WARNING: baseline model estimation failed")
        return(NULL)

    } else if( !inherits(fit.indep, "lavaan") ) {
        warning("lavaan WARNING: (user-provided) baseline model ",
                "is not a fitted lavaan object")
        return(NULL)

    } else if( !fit.indep@optim$converged ) {
        warning("lavaan WARNING: baseline model did not converge")
        return(NULL)

    } else {

        # evaluate if estimator/test matches original object
        # note: we do not need to check for 'se', as it may be 'none'
        sameTest <- all(object@Options$test == fit.indep@Options$test)
        if(!sameTest) {
            warning("lavaan WARNING:\n",
                    "\t Baseline model was using test(s) = ",
                    dQuote(fit.indep@Options$test),
                    "\n\t But original model was using test(s) = ",
                    dQuote(object@Options$test),
                    "\n\t Refitting baseline model!")
        }
        sameEstimator <- ( object@Options$estimator ==
                           fit.indep@Options$estimator )
        if(!sameEstimator) {
            warning("lavaan WARNING:\n",
                    "\t Baseline model was using estimator = ",
                    dQuote(fit.indep@Options$estimator),
                    "\n\t But original model was using estimator = ",
                    dQuote(object@Options$estimator),
                    "\n\t Refitting baseline model!")
        }
        if( !sameTest || !sameEstimator ) {
            lavoptions <- object@Options
            lavoptions$estimator   <- object@Options$estimator
            lavoptions$se          <- "none"
            lavoptions$verbose     <- FALSE
            lavoptions$baseline    <- FALSE
            lavoptions$check.start <- FALSE
            lavoptions$check.post  <- FALSE
            lavoptions$check.vcov  <- FALSE
            lavoptions$test        <- object@Options$test
            fit.indep <- try(lavaan(fit.indep,
                                    slotOptions     = lavoptions,
                                    slotData        = object@Data,
                                    slotSampleStats = object@SampleStats,
                                    sloth1          = object@h1,
                                    slotCache       = object@Cache),
                             silent = TRUE)
            # try again
            TEST <- lav_fit_measures_check_baseline(fit.indep = fit.indep,
                                                    object    = object)

        } else {
            # extract what we need
            TEST <- fit.indep@test
        }

    } # converged lavaan object

    TEST
}