R/lav_fit_srmr.R

In lavaan: Latent Variable Analysis

# functions related to the SRMR fit measures (single level only)

# lower-level functions:
# - lav_fit_srmr_mplus
# - lav_fit_srmr_twolevel

# higher-level functions:
# - lav_fit_srmr_lavobject

# Y.R. 22 July 2022

# Note: for rmrm/srmr/crmr, we use lav_residuals_summmary()

# SRMR for continuous data only
# see https://www.statmodel.com/download/SRMR.pdf
lav_fit_srmr_mplus <- function(lavobject) {

    # ngroups
    G <- lavobject@Data@ngroups

    # container per group
    srmr_mplus.group        <- numeric(G)
    srmr_mplus_nomean.group <- numeric(G)

    for(g in 1:G) {
        # observed
        if(!lavobject@SampleStats@missing.flag) {
            if(lavobject@Model@conditional.x) {
                S <- lavobject@SampleStats@res.cov[[g]]
                M <- lavobject@SampleStats@res.int[[g]]
            } else {
                S <- lavobject@SampleStats@cov[[g]]
                M <- lavobject@SampleStats@mean[[g]]
            }
        } else {
            # EM estimates
            S <- lavobject@SampleStats@missing.h1[[g]]$sigma
            M <- lavobject@SampleStats@missing.h1[[g]]$mu
        }
        nvar <- ncol(S)

        # estimated
        implied <- lavobject@implied
        lavmodel <- lavobject@Model
        Sigma.hat <- if(lavmodel@conditional.x) {
            implied$res.cov[[g]]
        } else {
            implied$cov[[g]]
        }
        Mu.hat <- if(lavmodel@conditional.x) {
            implied$res.int[[g]]
        } else {
            implied$mean[[g]]
        }

        # Bollen approach: simply using cov2cor ('residual correlations')
        S.cor <- cov2cor(S)
        Sigma.cor <- cov2cor(Sigma.hat)
        R.cor <- (S.cor - Sigma.cor)

        # meanstructure
        if(lavobject@Model@meanstructure) {
            # standardized residual mean vector
            R.cor.mean <- M/sqrt(diag(S)) - Mu.hat/sqrt(diag(Sigma.hat))

            e <- nvar*(nvar+1)/2 + nvar
            srmr_mplus.group[g] <-
                sqrt( (sum(R.cor[lower.tri(R.cor, diag = FALSE)]^2)  +
                       sum(R.cor.mean^2) +
                       sum(((diag(S) - diag(Sigma.hat))/diag(S))^2)) / e )

            e <- nvar*(nvar+1)/2
            srmr_mplus_nomean.group[g] <-
                sqrt( (sum(R.cor[lower.tri(R.cor, diag = FALSE)]^2)  +
                       sum(((diag(S) - diag(Sigma.hat))/diag(S))^2)) / e )
        } else {
            e <- nvar*(nvar+1)/2
            srmr_mplus_nomean.group[g] <- srmr_mplus.group[g] <-
                sqrt( (sum(R.cor[lower.tri(R.cor, diag = FALSE)]^2)  +
                       sum(((diag(S) - diag(Sigma.hat))/diag(S))^2)) / e )
        }
    } # G

    attr(srmr_mplus.group, "nomean") <- srmr_mplus_nomean.group
    srmr_mplus.group
}

lav_fit_srmr_twolevel <- function(lavobject = NULL) {

    nlevels <- lavobject@Data@nlevels
    G       <- lavobject@Data@ngroups

    SRMR.within  <- numeric(G)
    SRMR.between <- numeric(G)
    for(g in 1:G) {

        b.within  <- (g - 1L) * nlevels + 1L
        b.between <- (g - 1L) * nlevels + 2L

        # observed
        S.within  <- lavobject@h1$implied$cov[[  b.within  ]]
        M.within  <- lavobject@h1$implied$mean[[ b.within  ]]
        S.between <- lavobject@h1$implied$cov[[  b.between ]]
        M.between <- lavobject@h1$implied$mean[[ b.between ]]

        # estimated
        implied <- lav_model_implied_cond2uncond(lavobject@implied)
        Sigma.within  <- implied$cov[[  b.within  ]]
        Mu.within     <- implied$mean[[ b.within  ]]
        Sigma.between <- implied$cov[[  b.between ]]
        Mu.between    <- implied$mean[[ b.between ]]

        # force pd for between
        #    S.between <- lav_matrix_symmetric_force_pd(S.between)
        Sigma.between <- lav_matrix_symmetric_force_pd(Sigma.between)

        # Bollen approach: simply using cov2cor ('residual correlations')
        S.within.cor  <- cov2cor(S.within)
        S.between.cor <- cov2cor(S.between)
        Sigma.within.cor <- cov2cor(Sigma.within)
        if(all(diag(Sigma.between) > 0)) {
            Sigma.between.cor <- cov2cor(Sigma.between)
        } else {
            Sigma.between.cor <- matrix(as.numeric(NA),
                                     nrow = nrow(Sigma.between),
                                     ncol = ncol(Sigma.between))
        }
        R.within.cor  <- (S.within.cor  - Sigma.within.cor)
        R.between.cor <- (S.between.cor - Sigma.between.cor)

        nvar.within   <- NCOL(S.within)
        nvar.between  <- NCOL(S.between)
        pstar.within  <- nvar.within*(nvar.within+1)/2
        pstar.between <- nvar.between*(nvar.between+1)/2

        # SRMR
        SRMR.within[g] <-  sqrt( sum(lav_matrix_vech(R.within.cor)^2) /
                                 pstar.within )
        SRMR.between[g] <- sqrt( sum(lav_matrix_vech(R.between.cor)^2) /
                                 pstar.between )
    }

    # adjust for group sizes
    ng     <- unlist(lavobject@SampleStats@nobs)
    ntotal <- lavobject@SampleStats@ntotal
    SRMR_WITHIN  <- sum(ng/ntotal * SRMR.within)
    SRMR_BETWEEN <- sum(ng/ntotal * SRMR.between)
    SRMR_TOTAL <- SRMR_WITHIN + SRMR_BETWEEN

    c(SRMR_TOTAL, SRMR_WITHIN, SRMR_BETWEEN)
}

lav_fit_srmr_lavobject <- function(lavobject = NULL, fit.measures = "rmsea") {

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

    # categorical?
    categorical <- lavobject@Model@categorical

    # supported fit measures in this function
    if(categorical) {
        fit.srmr <- c("srmr")
        fit.srmr2 <- c("rmr", "rmr_nomean",
                       "srmr", # per default equal to srmr_bentler_nomean
                       "srmr_bentler", "srmr_bentler_nomean",
                       "crmr", "crmr_nomean",
                       "srmr_mplus", "srmr_mplus_nomean")
    } else {
        if(lavobject@Data@nlevels > 1L) {
            fit.srmr  <- c("srmr","srmr_within", "srmr_between")
            fit.srmr2 <- c("srmr","srmr_within", "srmr_between")
        } else {
            fit.srmr <- c("srmr")
            fit.srmr2 <- c("rmr", "rmr_nomean",
                           "srmr", # the default
                           "srmr_bentler", "srmr_bentler_nomean",
                           "crmr", "crmr_nomean",
                           "srmr_mplus", "srmr_mplus_nomean")
        }
    }

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

    # output container
    indices <- list()

    # 1. single level
    if(lavobject@Data@nlevels == 1L) {

        # RMR/SRMR/CRMR: we get it from lav_residuals_summary()
        out <- lav_residuals_summary(lavobject, se = FALSE, unbiased = FALSE)

        cov.cor <- "cov"
        if(categorical) {
            cov.cor <- "cor"
        }

        # only cov
         rmr_nomean.group  <- sapply(lapply(out, "[[",  "rmr"), "[[", cov.cor)
        srmr_nomean.group  <- sapply(lapply(out, "[[", "srmr"), "[[", cov.cor)
        crmr_nomean.group  <- sapply(lapply(out, "[[", "crmr"), "[[", cov.cor)

        # total
        if(lavobject@Model@meanstructure) {
             rmr.group <- sapply(lapply(out, "[[",  "rmr"), "[[", "total")
            srmr.group <- sapply(lapply(out, "[[", "srmr"), "[[", "total")
            crmr.group <- sapply(lapply(out, "[[", "crmr"), "[[", "total")
        } else {
            # no 'total', only 'cov'
             rmr.group <-  rmr_nomean.group
            srmr.group <- srmr_nomean.group
            crmr.group <- crmr_nomean.group
        }

        # the Mplus versions
        srmr_mplus.group <- lav_fit_srmr_mplus(lavobject = lavobject)
        srmr_mplus_nomean.group <- attr(srmr_mplus.group, "nomean")
        attr(srmr_mplus.group, "nomean") <- NULL

        # adjust for group sizes
        ng     <- unlist(lavobject@SampleStats@nobs)
        ntotal <- lavobject@SampleStats@ntotal
        RMR                 <- sum(ng/ntotal * rmr.group)
        RMR_NOMEAN          <- sum(ng/ntotal * rmr_nomean.group)
        SRMR_BENTLER        <- sum(ng/ntotal * srmr.group)
        SRMR_BENTLER_NOMEAN <- sum(ng/ntotal * srmr_nomean.group)
        CRMR                <- sum(ng/ntotal * crmr.group)
        CRMR_NOMEAN         <- sum(ng/ntotal * crmr_nomean.group)
        SRMR_MPLUS          <- sum(ng/ntotal * srmr_mplus.group)
        SRMR_MPLUS_NOMEAN   <- sum(ng/ntotal * srmr_mplus_nomean.group)

        # srmr
        if(lavobject@Options$mimic %in% c("lavaan", "EQS")) {
            if(categorical) {
                indices["srmr"] <- SRMR_BENTLER_NOMEAN
            } else {
                indices["srmr"] <- SRMR_BENTLER
            }
        } else if(lavobject@Options$mimic == "Mplus") {
            if(lavobject@Options$information[1] == "expected") {
                if(categorical) {
                    indices["srmr"] <- SRMR_BENTLER_NOMEAN
                } else {
                    indices["srmr"] <- SRMR_BENTLER
                }
            } else {
                if(categorical) {
                    indices["srmr"] <- SRMR_MPLUS_NOMEAN
                } else {
                   indices["srmr"] <- SRMR_MPLUS
                }
            }
        } # Mplus only

        # the others
        indices["srmr_bentler"]        <- SRMR_BENTLER
        indices["srmr_bentler_nomean"] <- SRMR_BENTLER_NOMEAN
        indices["crmr"]                <- CRMR
        indices["crmr_nomean"]         <- CRMR_NOMEAN

        # only correct for non-categorical:
        if(lavobject@Model@categorical) {
            # FIXME! Compute Mplus 8.1 way to compute SRMR in the
            #        categorical setting
            #        See 'SRMR in Mplus (2018)' document on Mplus website
            indices["srmr_mplus"]          <- as.numeric(NA)
            indices["srmr_mplus_nomean"]   <- as.numeric(NA)
        } else {
            indices["srmr_mplus"]          <- SRMR_MPLUS
            indices["srmr_mplus_nomean"]   <- SRMR_MPLUS_NOMEAN
        }
        indices["rmr"]             <- RMR
        indices["rmr_nomean"]          <- RMR_NOMEAN
    } else {

        # 2. twolevel setting
        out <- lav_fit_srmr_twolevel(lavobject = lavobject)
        indices["srmr"]         <- out[1]
        indices["srmr_within"]  <- out[2]
        indices["srmr_between"] <- out[3]

    } # twolevel

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