R/lav_residuals.R

In lavaan: Latent Variable Analysis

# residual diagnostics

# two types:
# 1) residuals for summary statistics
# 2) case-wise residuals

# this (new) version written around Aug/Sept 2018 for 0.6-3
# - based on obsList (inspect_sampstat) and estList (inspect_implied)
# - pre-scaling for type = "cor.bollen" and type = "cor.bentler"
# - summary statistics: rmr, srmr, crmr, urmr, usrmr, ucrmr; standard errors,
#       confidence intervals (for u(cs)rmr),
#       z-statistics (exact test, close test), p-values
# - type = "normalized" is based on lav_model_h1_acov(), and should now work
#   for all estimators
# - type = "standardized" now uses the correct formula, and should work for
#   for all estimators
# - type = "standardized.mplus" uses the simplified Mplus/LISREL version,
#   often resulting in NAs due to negative var(resid) estimates
#   (this was "standardized" in lavaan < 0.6.3

# WARNING: only partial support for conditional.x = TRUE!!
# - in categorical case: we only compute summary statistics, using cor + th
#   (no var, slopes, ...)
# - twolevel not supported here; see lav_fit_srmr.R, where we convert to
#   the unconditional setting

# - change 0.6-6: we enforce observed.information = "h1" to ensure 'Q' is a
#                 projection matrix (see lav_residuals_acov)

# - change 0.6-13: fixed.x = TRUE is ignored (to conform with 'tradition')

setMethod("residuals", "lavaan",
function(object, type = "raw", labels = TRUE) {

    # lowercase type
    type <- tolower(type)

    # type = "casewise"
    if(type %in% c("casewise","case","obs","observations","ov")) {
        return( lav_residuals_casewise(object, labels = labels) )
    } else {
        out <- lav_residuals(object = object, type = type, h1 = TRUE,
                             add.type = TRUE,
                             rename.cov.cor = FALSE, # should become FALSE!
                                                    # after packages (eg jmv)
                                                    # have adapted 0.6-3 style
                             add.labels = labels, add.class = TRUE,
                             drop.list.single.group = TRUE)
    }

    out
})

setMethod("resid", "lavaan",
function(object, type = "raw") {
    residuals(object, type = type, labels = TRUE)
})


# user-visible function
lavResiduals <- function(object, type = "cor.bentler", custom.rmr = NULL,
                         se = FALSE, zstat = TRUE, summary = TRUE,
                         h1.acov = "unstructured",
                         add.type = TRUE, add.labels = TRUE, add.class = TRUE,
                         drop.list.single.group = TRUE,
                         maximum.number = length(res.vech),
                         output = "list") {

    out <- lav_residuals(object = object, type = type, h1 = TRUE,
               custom.rmr = custom.rmr, se = se, zstat = zstat,
               summary = summary,
               summary.options = list(se = TRUE, zstat = TRUE, pvalue = TRUE,
                   unbiased = TRUE, unbiased.se = TRUE, unbiased.ci = TRUE,
                   unbiased.ci.level = 0.90, unbiased.zstat = TRUE,
                   unbiased.test.val = 0.05, unbiased.pvalue = TRUE),
               h1.acov = h1.acov, add.type = add.type,
               add.labels = add.labels, add.class = add.class,
               drop.list.single.group = drop.list.single.group)

    # no pretty printing yet...
    if(output == "table") {
        res <- out$cov
        # extract only below-diagonal elements
        res.vech <- lav_matrix_vech(res, diagonal = FALSE)

        # get names
        P <- nrow(res)
        NAMES <- colnames(res)

        nam <- expand.grid(NAMES,
                           NAMES)[lav_matrix_vech_idx(P, diagonal = FALSE),]
        NAMES.vech <- paste(nam[,1], "~~", nam[,2], sep = "")

        # create table
        TAB <- data.frame(name = NAMES.vech, res = round(res.vech, 3),
                          stringsAsFactors = FALSE)

        # sort table
        idx <- sort.int(abs(TAB$res), decreasing = TRUE,
                        index.return = TRUE)$ix
        out.sorted <- TAB[idx,]

        # show first rows only
        if(maximum.number == 0L || maximum.number > length(res.vech)) {
            maximum.number <- length(res.vech)
        }
        out <- out.sorted[seq_len(maximum.number), ]
    } else {
       # list -> nothing to do
    }

    out
}

# main function
lav_residuals <- function(object, type = "raw", h1 = TRUE, custom.rmr = NULL,
                          se = FALSE, zstat = FALSE, summary = FALSE,
                          summary.options = list(se = TRUE, zstat = TRUE,
                            pvalue = TRUE, unbiased = TRUE, unbiased.se = TRUE,
                            unbiased.ci = TRUE, unbiased.ci.level = 0.90,
                            unbiased.zstat = FALSE, unbiased.test.val = 0.05,
                            unbiased.pvalue = FALSE),
                          h1.acov = "unstructured", add.type = FALSE,
                          rename.cov.cor = FALSE,
                          add.labels = FALSE, add.class = FALSE,
                          drop.list.single.group = FALSE) {

    # type
    type <- tolower(type)[1]

    # check type
    if(!type %in% c("raw", "cor", "cor.bollen", "cor.bentler", "cor.eqs",
                    "rmr", "srmr", "crmr",
                    "normalized", "standardized", "standardized.mplus")) {
        stop("lavaan ERROR: unknown argument for type: ", dQuote(type))
    }

    # if cor, choose 'default'
    if(type == "cor") {
        if(object@Options$mimic == "EQS") {
            type <- "cor.bentler"
        } else {
            type <- "cor.bollen"
        }
    }
    if(type == "cor.eqs") {
        type <- "cor.bentler"
    }
    if(type == "rmr") {
        type <- "raw"
    }
    if(type == "srmr") {
        type <- "cor.bentler"
    }
    if(type == "crmr") {
        type <- "cor.bollen"
    }

    # slots
    lavdata  <- object@Data
    lavmodel <- object@Model

    # change options if multilevel (for now)
    if(lavdata@nlevels > 1L) {
        zstat <- se <- FALSE
        summary <- FALSE
    }

    # change options if categorical (for now)
    if(lavmodel@categorical) {

        # only if conditional.x = FALSE AND no continuous endogenous variables
        # -> only the simple setting where we only have thresholds and
        #    correlations

        # As soon as we add continuous variables, we get means/variances too,
        # and we need to decide how WLS.obs/WLS.est/WLS.V will then map to
        # the output of lavInspect(fit, "implied") and
        # lavInspect(fit, "sampstat")

        if(lavmodel@conditional.x || length(unlist(lavmodel@num.idx)) > 0L) {
            zstat <- se <- FALSE
            summary <- FALSE
            summary.options <- list(se = FALSE, zstat = FALSE,
                           pvalue = FALSE, unbiased = FALSE,
                           unbiased.se = FALSE,
                           unbiased.ci = FALSE, unbiased.ci.level = 0.90,
                           unbiased.zstat = FALSE, unbiased.test.val = 0.05,
                           unbiased.pvalue = FALSE)
        }
    }

    # change options if conditional.x (for now)
    if(!lavmodel@categorical && lavmodel@conditional.x) {
         zstat <- se <- FALSE
         summary <- FALSE
         summary.options <- list(se = FALSE, zstat = FALSE,
                            pvalue = FALSE, unbiased = FALSE,
                            unbiased.se = FALSE,
                            unbiased.ci = FALSE, unbiased.ci.level = 0.90,
                            unbiased.zstat = FALSE, unbiased.test.val = 0.05,
                            unbiased.pvalue = FALSE)
    }

    # observed and fitted sample statistics
    obsList <- lav_object_inspect_sampstat(object, h1 = h1,
                   add.labels = add.labels, add.class  = add.class,
                   drop.list.single.group = FALSE)
    estList <- lav_object_inspect_implied(object,
                   add.labels = add.labels, add.class = add.class,
                   drop.list.single.group = FALSE)
    # blocks
    nblocks <- length(obsList)

    # pre-scale?
    if(type %in% c("cor.bentler", "cor.bollen")) {
        for(b in seq_len(nblocks)) {
            var.obs <- if(lavmodel@conditional.x) {
                           diag(obsList[[b]][["res.cov"]])
                       } else {
                           diag(obsList[[b]][["cov"]])
                       }
            var.est <- if(lavmodel@conditional.x) {
                           diag(estList[[b]][["res.cov"]])
                       } else {
                           diag(estList[[b]][["cov"]])
                       }

            # rescale obsList
            obsList[[b]] <-
                lav_residuals_rescale(x = obsList[[b]], diag.cov = var.obs)
            # rescale estList
            if(type == "cor.bentler") { # use obsList
                estList[[b]] <-
                    lav_residuals_rescale(x = estList[[b]], diag.cov = var.obs)
            } else if(type == "cor.bollen") { # use estList for COV only
                estList[[b]] <- lav_residuals_rescale(x = estList[[b]],
                                    diag.cov  = var.est, diag.cov2 = var.obs)
            }
        }
    }

    # compute residuals: (observed - implied)
    resList <- vector("list", length = nblocks)
    for(b in seq_len(nblocks)) {
        resList[[b]] <- lapply(seq_len(length(obsList[[b]])),
                               FUN = function(el) {
                               obsList[[b]][[el]] - estList[[b]][[el]]})
        # always name the elements, even if add.labels = FALSE
        NAMES <- names(obsList[[b]])
        names(resList[[b]]) <- NAMES
    }

    # do we need seList?
    if(se || zstat) {
        seList <- lav_residuals_se(object, type = type, z.type = "standardized",
                               h1.acov = h1.acov,
                               add.class = add.class, add.labels = add.labels)
    } else if(type %in% c("normalized", "standardized", "standardized.mplus")) {
        seList <- lav_residuals_se(object, type = "raw", z.type = type,
                               h1.acov = h1.acov,
                               add.class = add.class, add.labels = add.labels)
    } else {
        seList <- NULL
    }

    # normalize/standardize?
    if(type %in% c("normalized", "standardized", "standardized.mplus")) {
        for(b in seq_len(nblocks)) {
            if(add.labels) {
                NAMES <- names(resList[[b]])
            }
            resList[[b]] <- lapply(seq_len(length(resList[[b]])),
                                   FUN = function(el) {
                                       A <- resList[[b]][[el]]
                                       B <-  seList[[b]][[el]]
                                       near.zero.idx <- which(abs(A) < 1e-05)
                                       if(length(near.zero.idx) > 0L) {
                                           B[near.zero.idx] <- 1
                                       }
                                       A/B })
            if(add.labels) {
                names(resList[[b]]) <- NAMES
            }
        }
    }

    # add se
    resList.orig <- resList
    if(se) {
        for(b in seq_len(nblocks)) {
            NAMES.res <- names(resList[[b]])
            NAMES.se  <- paste0(NAMES.res, ".se")
            resList[[b]] <- c(resList[[b]], seList[[b]])
            names(resList[[b]]) <- c(NAMES.res, NAMES.se)
        }
    }


    # add zstat
    if(zstat) {
        for(b in seq_len(nblocks)) {
            NAMES.res <- names(resList[[b]])
            NAMES.z   <- paste0(names(resList.orig[[b]]), ".z")
            tmp <- lapply(seq_len(length(resList.orig[[b]])),
                                   FUN = function(el) {
                                       A <- resList.orig[[b]][[el]]
                                       B <-  seList[[b]][[el]]
                                       # NOTE: which threshold should we use?
                                       # used to be 1e-05
                                       # changed to 1e-04 in 0.6-4
                                       near.zero.idx <- which(abs(A) < 1e-04)
                                       if(length(near.zero.idx) > 0L) {
                                           #B[near.zero.idx] <- as.numeric(NA)
                                           B[near.zero.idx] <- 1.0
                                       }
                                       A/B })
            resList[[b]] <- c(resList[[b]], tmp)
            names(resList[[b]]) <- c(NAMES.res, NAMES.z)
        }
    }

    # add summary statistics (rms, mabs)
    if(summary) {
        args <- c(list(object = object, type = type, h1.acov = h1.acov,
                       add.class = add.class, custom.rmr = custom.rmr),
                  summary.options)
        sumStat <- do.call("lav_residuals_summary", args)
        for(b in seq_len(nblocks)) {
            NAMES <- names(resList[[b]])
            resList[[b]] <- c(resList[[b]], list(sumStat[[b]][[1]])) # only 1
            NAMES <- c(NAMES, "summary")
            names(resList[[b]]) <- NAMES
        }
    }

    # last: add type
    if(add.type) {
        for(b in seq_len(nblocks)) {
            NAMES <- names(resList[[b]])
            resList[[b]] <- c(type, resList[[b]])
            NAMES <- c("type", NAMES)
            names(resList[[b]]) <- NAMES
        }
    }

    # optional: rename 'cov' to 'cor' (if type = "cor")
    if(rename.cov.cor && type %in% c("cor.bentler", "cor.bollen")) {
        for(b in seq_len(nblocks)) {
            NAMES <- names(resList[[b]])
            NAMES <- gsub("cov", "cor", NAMES)
            names(resList[[b]]) <- NAMES
        }
    }


    # output
    OUT <- resList
    if(nblocks == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(lavdata@nlevels == 1L &&
           length(lavdata@group.label) > 0L) {
            names(OUT) <- unlist(lavdata@group.label)
        } else if(lavdata@nlevels > 1L &&
                  length(lavdata@group.label) == 0L) {
            names(OUT) <- lavdata@level.label
        }
    }

    OUT
}

# return ACOV as list per group
lav_residuals_acov <- function(object, type = "raw", z.type = "standardized",
                               h1.acov = "unstructured") {

    # check type
    if(z.type %in% c("normalized", "standardized.mplus") && type != "raw") {
        stop("lavaan ERROR: z.type = ", dQuote(z.type), " can only be used ",
             "with type = ", dQuote("raw"))
    }

    # slots
    lavdata        <- object@Data
    lavmodel       <- object@Model
    lavsamplestats <- object@SampleStats

    # return list per group
    ACOV.res <- vector("list", length = lavdata@ngroups)

    # compute ACOV for observed h1 sample statistics (ACOV == Gamma/N)
    if(lavmodel@categorical) {
        NACOV.obs <- lavsamplestats@NACOV
        ACOV.obs <- lapply(NACOV.obs, function(x) x / lavsamplestats@ntotal)
    } else {
        ACOV.obs <- lav_model_h1_acov(lavobject = object,
                                      h1.information = h1.acov)
    }

    # shortcut for normalized
    if(z.type == "normalized") {
        ACOV.res <- ACOV.obs
        return(ACOV.res)
    } else {
        if(z.type == "standardized") {
            A1 <- lav_model_h1_information(object)
            if(lavmodel@estimator == "DWLS" || lavmodel@estimator == "ULS") {
               # A1 is diagonal matrix
               A1 <- lapply(A1, diag)
            }
            if(type %in% c("cor.bentler", "cor.bollen")) {
                sampstat <- lavTech(object, "sampstat")
            }
        } else if(z.type == "standardized.mplus") {
            VCOV  <- lavTech(object, "vcov")
        }
        DELTA <- lavTech(object, "delta")
    }

    # for each group, compute ACOV
    for(g in seq_len(lavdata@ngroups)) {

        # group weight
        gw <- object@SampleStats@nobs[[g]] / object@SampleStats@ntotal

        if(z.type == "standardized.mplus") { # simplified formula
                                             # also used by LISREL?
            # see https://www.statmodel.com/download/StandardizedResiduals.pdf

            ACOV.est.g <- DELTA[[g]] %*% VCOV %*% t(DELTA[[g]])
            ACOV.res[[g]] <- ACOV.obs[[g]] - ACOV.est.g

        } else if(z.type == "standardized") {
            # see Ogasawara (2001) using Bentler & Dijkstra (1985) eq 1.7.4

            # NVarCov, but always 'not' robust
            #
            # new in 0.6-6: to ensure Q is a projection matrix, we
            #               force observed.information = "h1"
            #               (only needed if information is observed)
            this.options <- object@Options
            this.options$observed.information[1] <- "h1"
            A0.g.inv <- lav_model_information(lavmodel  = lavmodel,
                                         lavsamplestats = object@SampleStats,
                                         lavdata        = lavdata,
                                         lavcache       = object@Cache,
                                         lavimplied     = object@implied,
                                         lavh1          = object@h1,
                                         lavoptions     = this.options,
                                         extra          = FALSE,
                                         augmented      = TRUE,
                                         inverted       = TRUE,
                                         use.ginv       = TRUE)

            ACOV.est.g <- gw * (DELTA[[g]] %*% A0.g.inv %*% t(DELTA[[g]]))
            Q <- diag(nrow = nrow(ACOV.est.g)) - ACOV.est.g %*% A1[[g]]
            ACOV.res[[g]] <- Q %*% ACOV.obs[[g]] %*% t(Q)

            # correct ACOV.res for type = "cor.bentler" or type = "cor.bollen"
            if(type == "cor.bentler") {
                if(lavmodel@categorical) {
                    if(lavmodel@conditional.x ||
                       length(unlist(lavmodel@num.idx)) > 0L) {
                    stop("lavaan ERROR: SE for cor.bentler not available (yet) if categorical = TRUE, and conditional.x = TRUE OR some endogenous variables are continuous")        } else {
                       # nothing to do, as we already are in correlation metric
                    }

                } else {
                    # Ogasawara (2001), eq (13), or
                    # Maydeu-Olivares (2017), eq (16)
                    COV <- if(lavmodel@conditional.x) {
                               sampstat[[g]][["res.cov"]]
                           } else { sampstat[[g]][["cov"]] }
                    SS <- 1/sqrt(diag(COV))
                    tmp <- lav_matrix_vech(tcrossprod(SS))
                    G.inv.sqrt <- diag( tmp, nrow = length(tmp) )
                    if(lavmodel@meanstructure) {
                        GG <- lav_matrix_bdiag(diag(SS, nrow = length(SS)),
                                               G.inv.sqrt)
                    } else {
                        GG <- G.inv.sqrt
                    }
                    ACOV.res[[g]] <- GG %*% ACOV.res[[g]] %*% GG
                } # continuous


            } else if(type == "cor.bollen") {
                if(lavmodel@categorical) {
                    if(lavmodel@conditional.x ||
                       length(unlist(lavmodel@num.idx)) > 0L) {
                    stop("lavaan ERROR: SE for cor.bentler not available (yet) if categorical = TRUE, and conditional.x = TRUE OR some endogenous variables are continuous")        } else {
                       # nothing to do, as we already are in correlation metric
                    }

                } else {
                    # here we use the Maydeu-Olivares (2017) approach, see eq 17
                    COV <- if(lavmodel@conditional.x) {
                               sampstat[[g]][["res.cov"]]
                           } else { sampstat[[g]][["cov"]] }
                    F1 <- lav_deriv_cov2corB(COV)
                    if(lavmodel@meanstructure) {
                        SS <- 1/sqrt(diag(COV))
                        FF <- lav_matrix_bdiag(diag(SS, nrow = length(SS)), F1)
                    } else {
                        FF <- F1
                    }
                    ACOV.res[[g]] <- FF %*% ACOV.res[[g]] %*% t(FF)
                } # continuous
            } # cor.bollen

        } # z.type = "standardized"

    } # g

    ACOV.res
}

# return resList with 'se' values for each residual
lav_residuals_se <- function(object, type = "raw", z.type = "standardized",
                             h1.acov = "unstructured",
                             add.class = FALSE, add.labels = FALSE) {

    # slots
    lavdata    <- object@Data
    lavmodel   <- object@Model
    lavpta     <- object@pta

    # return list per group
    seList <- vector("list", length = lavdata@ngroups)

    # get ACOV per group
    ACOV.res <- lav_residuals_acov(object = object, type = type,
                                   z.type = z.type, h1.acov = h1.acov)

    # labels
    if(add.labels) {
        ov.names <- object@pta$vnames$ov
        ov.names.res <- object@pta$vnames$ov.nox
        ov.names.x   <- object@pta$vnames$ov.x
    }

    # for each group, compute 'se' values, and fill list
    for(g in seq_len(lavdata@ngroups)) {

        nvar <- object@pta$nvar[[g]] # block or group-based?
        diag.ACOV <- diag(ACOV.res[[g]])

        # take care of negative, or non-finite diag.ACOV elements
        diag.ACOV[!is.finite(diag.ACOV)] <- NA
        diag.ACOV[ diag.ACOV < 0 ] <- NA

        # categorical
        if(lavmodel@categorical) {
            if(lavmodel@conditional.x ||
               length(unlist(lavmodel@num.idx)) > 0L) {
                stop("not ready yet!")
            }

            # COR
            nth <- length(lavmodel@th.idx[[g]])
            tmp <- sqrt(diag.ACOV[-(1:nth)])
            cov.se <- lav_matrix_vech_reverse(tmp, diagonal = FALSE)

            # MEAN
            mean.se <- rep(as.numeric(NA), nth)

            # TH
            th.se <- sqrt(diag.ACOV[1:nth])

            if(add.class) {
                class(cov.se) <- c("lavaan.matrix.symmetric", "matrix")
                class(mean.se) <- c("lavaan.vector", "numeric")
                class(th.se) <- c("lavaan.vector", "numeric")
            }
            if(add.labels) {
                rownames(cov.se) <- colnames(cov.se) <- ov.names[[g]]
                names(mean.se) <- ov.names[[g]]
                names(th.se) <- lavpta$vnames$th.mean[[g]]
            }
            seList[[g]] <- list(cov.se = cov.se, mean.se = mean.se,
                                th.se = th.se)


        # continuous -- single level
        } else if(lavdata@nlevels == 1L) {
            if(lavmodel@conditional.x) {
                stop("not ready yet")
            } else {
                if(lavmodel@meanstructure) {
                    tmp <- sqrt(diag.ACOV[-(1:nvar)])
                    cov.se <- lav_matrix_vech_reverse(tmp)
                    mean.se <- sqrt(diag.ACOV[1:nvar])
                    if(add.class) {
                        class(cov.se) <- c("lavaan.matrix.symmetric", "matrix")
                        class(mean.se) <- c("lavaan.vector", "numeric")
                    }
                    if(add.labels) {
                        rownames(cov.se) <- colnames(cov.se) <- ov.names[[g]]
                        names(mean.se) <- ov.names[[g]]
                    }
                    seList[[g]] <- list(cov.se = cov.se, mean.se = mean.se)
                } else {
                    cov.se <- lav_matrix_vech_reverse(sqrt(diag.ACOV))
                    if(add.class) {
                        class(cov.se) <- c("lavaan.matrix.symmetric", "matrix")
                    }
                    if(add.labels) {
                        rownames(cov.se) <- colnames(cov.se) <- ov.names[[g]]
                    }
                    seList[[g]] <- list(cov.se = cov.se)
                }
            }

        # continuous -- multilevel
        } else if(lavdata@nlevels > 1L) {
            stop("not ready yet")
        }

    } # g

    seList
}

# return summary statistics as list per group
lav_residuals_summary <- function(object, type = c("rmr", "srmr", "crmr"),
                                  h1.acov = "unstructured", custom.rmr = NULL,
                                  se = FALSE, zstat = FALSE, pvalue = FALSE,
                                  unbiased = FALSE, unbiased.se = FALSE,
                                  unbiased.ci = FALSE, unbiased.ci.level = 0.90,
                                  unbiased.zstat = FALSE,
                                  unbiased.test.val = 0.05,
                                  unbiased.pvalue = FALSE,
                                  add.class = FALSE) {

    # arguments
    if (length(custom.rmr)) {
        if (!is.list(custom.rmr)) stop('custom.rmr must be a list')
        ## Each custom (S/C)RMR must have a unique name
        customNAMES <- names(custom.rmr)
        if (is.null(customNAMES)) stop('custom.rmr list must have names')
        if (length(unique(customNAMES)) < length(custom.rmr)) {
            stop('custom.rmr must have a unique name for each summary')
        }
        ## Each list must contain a list consisting of $cov and/or $mean (no $th yet)
        for (i in seq_along(custom.rmr)) {
            if (!is.list(custom.rmr[[i]])) {
                stop('Each element in custom.rmr must be a list')
            }
            if (is.null(names(custom.rmr[[i]]))) {
                stop('The list in custom.rmr must have names')
            }
            if (!all(names(custom.rmr[[i]]) %in% c("cov","mean"))) {
                stop('Elements in custom.rmr must be names "cov" and/or "mean"')
            }
            ## below, verify dimensions match rmsList.g
        }
        #FIXME: blocks can have unique models, need another layer of lists
        #       between custom summaries and moments
    } else {
        customNAMES <- NULL
    }

    if(pvalue) {
        zstat <- TRUE
    }
    if(zstat) {
        se <- TRUE
    }
    if(unbiased.pvalue) {
        unbiased.zstat <- TRUE
    }
    if(unbiased.zstat) {
        unbiased.se <- TRUE
    }

    if(!all(type %in% c("rmr", "srmr", "crmr",
                        "raw", "cor.bentler", "cor.bollen"))) {
        stop("lavaan ERROR: unknown type: ", dQuote(type))
    }

    # change type name
    idx <- which(type == "raw")
    if(length(idx) > 0L) {
        type[idx] <- "rmr"
    }
    idx <- which(type == "cor.bentler")
    if(length(idx) > 0L) {
        type[idx] <- "srmr"
    }
    idx <- which(type == "cor.bollen")
    if(length(idx) > 0L) {
        type[idx] <- "crmr"
    }

    # slots
    lavdata        <- object@Data
    lavmodel       <- object@Model

    # fixed.x/conditional.x
    fixed.x <- lavmodel@fixed.x
    conditional.x <- lavmodel@conditional.x


    rmrFlag <- srmrFlag <- crmrFlag <- FALSE
    if("rmr" %in% type || "raw" %in% type) {
        # FIXME: recursive call to lav_residuals() is summary = TRUE!!
        rmrList <- lav_residuals(object = object, type = "raw")
        if(se || unbiased) {
            rmrList.se <- lav_residuals_acov(object = object, type = "raw",
                                             z.type = "standardized",
                                             h1.acov = "unstructured")
        }
    }
    if("srmr" %in% type || "cor.bentler" %in% type || "cor" %in% type) {
        srmrList <- lav_residuals(object = object, type = "cor.bentler")
        if(se || unbiased) {
            srmrList.se <- lav_residuals_acov(object = object,
                                              type = "cor.bentler",
                                              z.type = "standardized",
                                              h1.acov = "unstructured")
        }
    }
    if("crmr" %in% type || "cor.bollen" %in% type) {
        crmrList <- lav_residuals(object = object, type = "cor.bollen")
        if(se || unbiased) {
            crmrList.se <- lav_residuals_acov(object = object,
                                              type = "cor.bollen",
                                              z.type = "standardized",
                                              h1.acov = "unstructured")
        }
    }

    # return list per group
    sumStat <- vector("list", length = lavdata@ngroups)

    # for each group, compute ACOV
    for(g in seq_len(lavdata@ngroups)) {

        nvar <- object@pta$nvar[[g]] # block or group-based?

        # categorical single level
        if(lavdata@nlevels == 1L && lavmodel@categorical) {

            if((se || unbiased) && (conditional.x ||
               length(unlist(lavmodel@num.idx)) > 0L)) {
                stop("not ready yet")
            } else {

                # remove fixed.x elements:
                # seems like a good idea, but nobody likes it
                # nvar.x <- pstar.x <- 0L
                # if(lavmodel@fixed.x) {
                #     nvar.x <- lavmodel@nexo[g]
                #     pstar.x <- nvar.x * (nvar.x - 1) / 2 # note '-'
                # }

                OUT <- vector("list", length(type))
                names(OUT) <- type

                for(typ in seq_len(length(type))) {

                    if(type[typ] == "rmr") {
                        rmsList.g <- rmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- rmrList.se[[g]]
                        }
                    } else if(type[typ] == "srmr") {
                        rmsList.g <- srmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- srmrList.se[[g]]
                        }
                    } else if(type[typ] == "crmr") {
                        rmsList.g <- crmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- crmrList.se[[g]]
                        }
                    }

                    # COR
                    nth <- length(lavmodel@th.idx[[g]])
                    if(conditional.x) {
                        STATS <- lav_matrix_vech(rmsList.g[["res.cov"]],
                                                 diagonal = FALSE)
                    } else {
                        STATS <- lav_matrix_vech(rmsList.g[["cov"]],
                                                 diagonal = FALSE)
                    }

                    # should pstar be p*(p+1)/2 or p*(p-1)/2
                    # we use the first for SRMR and the latter for CRMR
                    if(type[typ] == "crmr") {
                        pstar <- length(STATS)
                    } else {
                        pstar <- length(STATS) + nvar
                    }
                    ACOV <- NULL
                    if(se || unbiased) {
                        ACOV <- rmsList.se.g[-seq_len(nth),
                                             -seq_len(nth), drop = FALSE]
                    }
                    RMS.COR <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])


                    # THRESHOLDS
                    if(conditional.x) {
                        STATS <- rmsList.g[["res.th"]]
                    } else {
                        STATS <- rmsList.g[["th"]]
                    }
                    pstar <- length(STATS)
                    ACOV <- NULL
                    if(se || unbiased) {
                        ACOV <- rmsList.se.g[seq_len(nth),
                                             seq_len(nth), drop = FALSE]
                    }
                    RMS.TH <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])

                    # MEAN
                    #STATS <- rmsList.g[["mean"]]
                    STATS <- numeric(0L)
                    pstar <- length(STATS)
                    ACOV <- NULL
                    if(se || unbiased) {
                        # TODO: extract from rmsList.se.g
                    }
                    RMS.MEAN <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])

                    # VAR (not ready yet)
                    #STATS <- diag(rmsList.g[["cov"]])[lavmodel@num.idx[[g]]]
                    STATS <- numeric(0L)
                    pstar <- length(STATS)
                    ACOV <- NULL
                    if(se || unbiased) {
                        # TODO: extract from rmsList.se.g
                    }
                    RMS.VAR <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])

                    # TOTAL -- FIXME: for conditional.x ....
                    if(conditional.x) {
                        STATS <- c(lav_matrix_vech(rmsList.g[["res.cov"]],
                                                   diagonal = FALSE),
                                   rmsList.g[["res.th"]])
                    } else {
                        STATS <- c(lav_matrix_vech(rmsList.g[["cov"]],
                                                   diagonal = FALSE),
                                   rmsList.g[["th"]])
                                   #rmsList.g[["mean"]],
                               #diag(rmsList.g[["cov"]])[lavmodel@num.idx[[g]]])
                    }

                    # should pstar be p*(p+1)/2 or p*(p-1)/2 for COV/COR?
                    # we use the first for SRMR and the latter for CRMR
                    if(type[typ] == "crmr") {
                        pstar <- length(STATS)
                    } else {
                        pstar <- length(STATS) + nvar
                    }

                    #if(lavmodel@fixed.x) {
                    #    pstar <- pstar - pstar.x
                    #}

                    ACOV <- NULL
                    if(se || unbiased) {
                        ACOV <- rmsList.se.g
                    }
                    RMS.TOTAL <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])

                    TABLE <- as.data.frame(cbind(RMS.COR,
                                                 RMS.TH,
                                                # RMS.MEAN,
                                                # RMS.VAR,
                                                 RMS.TOTAL))
                    #colnames(TABLE) <- c("cor", "thresholds", "mean",
                    #                     "var", "total")
                    colnames(TABLE) <- c("cor", "thresholds", "total")
                    if(add.class) {
                        class(TABLE) <- c("lavaan.data.frame", "data.frame")
                    }
                    OUT[[typ]] <- TABLE

                } # type

            } # not conditional.x or mixed cat/con

        # continuous -- single level
        } else if(lavdata@nlevels == 1L) {
            if((se || unbiased) && conditional.x) {
                stop("not ready yet")
            } else {

                #nvar.x <- pstar.x <- 0L
                #if(lavmodel@fixed.x) {
                #    nvar.x <- lavmodel@nexo[g]
                #    pstar.x <- nvar.x * (nvar.x + 1) / 2
                #}

                OUT <- vector("list", length(type))
                names(OUT) <- type

                for(typ in seq_len(length(type))) {

                    if(type[typ] == "rmr") {
                        rmsList.g <- rmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- rmrList.se[[g]]
                        }
                    } else if(type[typ] == "srmr") {
                        rmsList.g <- srmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- srmrList.se[[g]]
                        }
                    } else if(type[typ] == "crmr") {
                        rmsList.g <- crmrList[[g]]
                        if(se || unbiased) {
                            rmsList.se.g <- crmrList.se[[g]]
                        }
                    }

                    # COV
                    if(conditional.x) {
                        STATS <- lav_matrix_vech(rmsList.g[["res.cov"]])
                    } else {
                        STATS <- lav_matrix_vech(rmsList.g[["cov"]])
                    }
                    #pstar <- ( length(STATS) - pstar.x )
                    pstar <- length(STATS)
                    if(type[typ] == "crmr") {
                        # pstar <- pstar - ( nvar - nvar.x )
                        pstar <- pstar - nvar
                    }

                    ACOV <- NULL
                    if(se || unbiased) {
                        ACOV <- if(lavmodel@meanstructure) {
                                  rmsList.se.g[-seq_len(nvar),
                                               -seq_len(nvar), drop = FALSE]
                                } else { rmsList.se.g }
                    }
                    RMS.COV <- lav_residuals_summary_rms(STATS = STATS,
                        ACOV = ACOV, se = se, zstat = zstat, pvalue = pvalue,
                        unbiased = unbiased, unbiased.se = unbiased.se,
                        unbiased.ci = unbiased.ci,
                        unbiased.ci.level = unbiased.ci.level,
                        unbiased.zstat = unbiased.zstat,
                        unbiased.test.val = unbiased.test.val,
                        unbiased.pvalue = unbiased.pvalue,
                        pstar = pstar, type = type[typ])

                    # MEAN
                    if(lavmodel@meanstructure) {
                        if(conditional.x) {
                            STATS <- rmsList.g[["res.int"]]
                        } else {
                            STATS <- rmsList.g[["mean"]]
                        }
                        # pstar <- ( length(STATS) - nvar.x )
                        pstar <- length(STATS)
                        ACOV  <- NULL
                        if(se || unbiased) {
                            ACOV <- rmsList.se.g[seq_len(nvar),
                                                 seq_len(nvar), drop = FALSE]
                        }
                        RMS.MEAN <- lav_residuals_summary_rms(STATS = STATS,
                            ACOV = ACOV,
                            se = se, zstat = zstat, pvalue = pvalue,
                            unbiased = unbiased, unbiased.se = unbiased.se,
                            unbiased.ci = unbiased.ci,
                            unbiased.ci.level = unbiased.ci.level,
                            unbiased.zstat = unbiased.zstat,
                            unbiased.test.val = unbiased.test.val,
                            unbiased.pvalue = unbiased.pvalue,
                            pstar = pstar, type = type[typ])
                    }

                    # TOTAL
                    if(lavmodel@meanstructure) {
                        if(conditional.x) {
                            STATS <- c(rmsList.g[["res.int"]],
                                       lav_matrix_vech(rmsList.g[["res.cov"]]))
                        } else {
                            STATS <- c(rmsList.g[["mean"]],
                                       lav_matrix_vech(rmsList.g[["cov"]]))
                        }
                        # pstar <- ( length(STATS) - ( pstar.x + nvar.x) )
                        pstar <- length(STATS)
                        if(type[typ] == "crmr") {
                            # pstar <- pstar - ( nvar - nvar.x )
                            pstar <- pstar - nvar
                        }
                        ACOV  <- NULL
                        if(se || unbiased) {
                            ACOV <- rmsList.se.g
                        }
                        RMS.TOTAL <- lav_residuals_summary_rms(STATS = STATS,
                            ACOV = ACOV,
                            se = se, zstat = zstat, pvalue = pvalue,
                            unbiased = unbiased, unbiased.se = unbiased.se,
                            unbiased.ci = unbiased.ci,
                            unbiased.ci.level = unbiased.ci.level,
                            unbiased.zstat = unbiased.zstat,
                            unbiased.test.val = unbiased.test.val,
                            unbiased.pvalue = unbiased.pvalue,
                            pstar = pstar, type = type[typ])
                    }

                    # CUSTOM
                    if (length(custom.rmr)) {
                      if (lavmodel@fixed.x && !lavmodel@conditional.x) {
                        ## save exogenous-variable indices, use to remove or set
                        ## FALSE any moments that cannot have nonzero residuals
                        x.idx <- which(rownames(rmsList.g$cov) %in% object@Data@ov.names.x[[g]])
                      }

                      RMS.CUSTOM.LIST <- vector("list", length(customNAMES))

                      for (cus in customNAMES) {
                        ## in case there is no meanstructure
                        STATS <- NULL
                        ACOV.idx  <- NULL

                        # MEANS?
                        if (lavmodel@meanstructure) {
                          if ("mean" %in% names(custom.rmr[[cus]])) {

                            ## if logical, save numeric indices
                            if (is.logical(custom.rmr[[cus]]$mean)) {
                              ## check length
                              if (length(custom.rmr[[cus]]$mean) != length(rmsList.g[["mean"]])) {
                                stop('length(custom.rmr$', cus, '$mean) must ',
                                     'match length(lavResiduals(fit)$mean)')
                              }
                              ACOV.idx <- which(custom.rmr[[cus]]$mean)
                              if (lavmodel@fixed.x && !lavmodel@conditional.x) {
                                ACOV.idx[x.idx] <- FALSE
                              }

                            } else if (!is.numeric(custom.rmr[[cus]]$mean)) {
                              stop('custom.rmr$', cus, '$mean must contain ',
                                   'logical or numeric indices.')

                            } else {
                              ACOV.idx <- custom.rmr[[cus]]$mean
                              if (lavmodel@fixed.x && !lavmodel@conditional.x) {
                                  ACOV.idx <- setdiff(ACOV.idx, x.idx)
                              }
                              ACOV.idx <- ACOV.idx[!is.na(ACOV.idx)] # necessary?
                              if (max(ACOV.idx) > length(rmsList.g[["mean"]])) {
                                stop('custom.rmr$', cus, '$mean[',
                                     which.max(ACOV.idx),
                                     '] is an out-of-bounds index')
                              }
                            }

                            STATS <- rmsList.g[["mean"]][ACOV.idx]
                          }
                        }


                        # (CO)VARIANCES?
                        if ("cov" %in% names(custom.rmr[[cus]])) {

                          ## if numeric, create a logical matrix to obtain
                          ## ACOV.idx and check for x.idx
                          if (is.numeric(custom.rmr[[cus]]$cov)) {
                            cusCOV <- rmsList.g[["cov"]] == "start with all FALSE"
                            ## matrix of row/column indices?
                            if (length(dim(custom.rmr[[cus]]$cov))) {
                              if (max(custom.rmr[[cus]]$cov[,1:2] > nrow(rmsList.g[["cov"]]))) {
                                stop('numeric indices in custom.rmr$', cus, '$cov',
                                     ' cannot exceed ', nrow(rmsList.g[["cov"]]))
                              }
                              for (RR in 1:nrow(custom.rmr[[cus]]$cov)) {
                                cusCOV[ custom.rmr[[cus]]$cov[RR, 1] ,
                                        custom.rmr[[cus]]$cov[RR, 2] ] <- TRUE
                              }

                            } else {
                              ## numeric-vector indices
                              if (max(custom.rmr[[cus]]$cov > length(rmsList.g[["cov"]]))) {
                                stop('numeric indices in custom.rmr$', cus, '$cov',
                                     ' cannot exceed ', length(rmsList.g[["cov"]]))
                              }
                              cusCOV[custom.rmr[[cus]]$cov] <- TRUE

                            }

                            ## numeric indices no longer needed, use logical
                            custom.rmr[[cus]]$cov <- cusCOV

                          } else if (!is.logical(custom.rmr[[cus]]$cov)) {
                              stop('custom.rmr$', cus, '$cov must be a logical ',
                                   'square matrix or a numeric matrix of ',
                                   '(row/column) indices.')
                          }

                          ## check dimensions
                          if (!all(dim(custom.rmr[[cus]]$cov) == dim(rmsList.g[["cov"]]))) {
                            stop('dim(custom.rmr$', cus, '$cov) must ',
                                 'match dim(lavResiduals(fit)$cov)')
                          }
                          ## users can specify upper.tri or lower.tri indices
                          custom.rmr[[cus]]$cov <- custom.rmr[[cus]]$cov | t(custom.rmr[[cus]]$cov)
                          ## but ACOV refers to lower.tri indices
                          custom.rmr[[cus]]$cov[upper.tri(custom.rmr[[cus]]$cov)] <- FALSE
                          ## diagonal relevant?
                          if (type[typ] == "crmr") diag(custom.rmr[[cus]]$cov) <- FALSE
                          ## extract lower.tri indices
                          vech.idx <- which(lav_matrix_vech(custom.rmr[[cus]]$cov))

                          ## add residuals to STATS, indices to ACOV.idx
                          STATS <- c(STATS, lav_matrix_vech(rmsList.g[["cov"]])[vech.idx])
                          ACOV.idx <- c(ACOV.idx, vech.idx)

                        }


                        ## count residuals in summary (x.idx already removed)
                        pstar <- length(STATS)

                        ACOV  <- NULL
                        if (se || unbiased) {
                          ACOV <- rmsList.se.g[ACOV.idx, ACOV.idx, drop = FALSE]
                        }
                        RMS.CUSTOM.LIST[[cus]] <- lav_residuals_summary_rms(STATS = STATS,
                            ACOV = ACOV,
                            se = se, zstat = zstat, pvalue = pvalue,
                            unbiased = unbiased, unbiased.se = unbiased.se,
                            unbiased.ci = unbiased.ci,
                            unbiased.ci.level = unbiased.ci.level,
                            unbiased.zstat = unbiased.zstat,
                            unbiased.test.val = unbiased.test.val,
                            unbiased.pvalue = unbiased.pvalue,
                            pstar = pstar, type = type[typ])

                      #FIXME: update for categorical
                      } # cus
                        RMS.CUSTOM <- do.call(rbind, RMS.CUSTOM.LIST)
                    } else {
                        RMS.CUSTOM <- NULL
                    }


                    if(lavmodel@meanstructure) {
                        TABLE <- as.data.frame(cbind(RMS.COV,
                                                     RMS.MEAN,
                                                     RMS.TOTAL,
                                                     RMS.CUSTOM))
                        colnames(TABLE) <- c("cov", "mean", "total",
                                             customNAMES)
                    } else {
                        TABLE <- as.data.frame(cbind(RMS.COV, RMS.CUSTOM))
                        colnames(TABLE) <- c("cov", customNAMES)
                    }
                    if(add.class) {
                        class(TABLE) <- c("lavaan.data.frame", "data.frame")
                    }
                    OUT[[typ]] <- TABLE

                } # type

            } # continuous, single-level, unconditional

        # continuous -- multilevel
        } else if(lavdata@nlevels > 1L) {
            stop("not ready yet")
        }

        sumStat[[g]] <- OUT
    } # g

    sumStat
}


lav_residuals_summary_rms <- function(STATS = NULL, ACOV = NULL,
                                      se = FALSE,
                                      level = 0.90,
                                      zstat = FALSE, pvalue = FALSE,
                                      unbiased = FALSE,
                                      unbiased.se = FALSE,
                                      unbiased.ci = FALSE,
                                      unbiased.ci.level = 0.90,
                                      unbiased.zstat = FALSE,
                                      unbiased.test.val = 0.05,
                                      unbiased.pvalue = FALSE,
                                      pstar = 0, type = "rms") {

    OUT <- vector("list", length = 0L)


    # covariance matrix
    if(length(STATS) > 0L) {
        rms <- sqrt(sum(STATS * STATS)/pstar)
    } else {
        rms <- 0
        se <- unbiased <- zstat <- FALSE
    }

    # default is NULL
    rms.se <- rms.z <- rms.pvalue <- NULL
    urms <- urms.se <- urms.z <- urms.pvalue <- NULL
    urms.ci.lower <- urms.ci.upper <- NULL
    if(!unbiased.zstat) {
        unbiased.test.val <- NULL
    }

    if(se || unbiased) {
        TR2 <- sum(diag( ACOV %*% ACOV ))
        TR1 <- sum(diag( ACOV ))
        if(se) {
            rms.avar <- TR2/(TR1 * 2 * pstar)
            if(!is.finite(rms.avar) || rms.avar < .Machine$double.eps) {
                rms.se <- as.numeric(NA)
            } else {
                rms.se <- sqrt(rms.avar)
            }
        }
    }

    if(zstat) {
        E.rms <- ( sqrt(TR1/pstar) * (4 * TR1 * TR1 - TR2) / (4 * TR1 * TR1) )
        rms.z <- max((rms - E.rms), 0) / rms.se
        if(pvalue) {
            rms.pvalue <- 1 - pnorm(rms.z)
        }
    }

    if(unbiased) {
        T.cov <- as.numeric( crossprod(STATS) )
        eVe <- as.numeric(t(STATS) %*% ACOV %*% STATS)
        k.cov <- 1 - (TR2 + 2 * eVe) / (4 * T.cov * T.cov)
        urms <- ( 1/k.cov * sqrt( max((T.cov - TR1), 0) / pstar ) )
        if(unbiased.se) {
            urms.avar <- ( 1/(k.cov*k.cov) * (TR2 + 2*eVe) / (2*pstar * T.cov) )
            if(!is.finite(urms.avar) || urms.avar < .Machine$double.eps) {
                urms.se <- as.numeric(NA)
            } else {
                urms.se <- sqrt(urms.avar)
            }
            if(unbiased.ci) {
                a <- (1 - unbiased.ci.level)/2
                a <- c(a, 1-a)
                fac <- stats::qnorm(a)
                urms.ci.lower <- urms + urms.se * fac[1]
                urms.ci.upper <- urms + urms.se * fac[2]
            }
            if(unbiased.zstat) {
                urms.z <- (urms - unbiased.test.val) / urms.se
                if(unbiased.pvalue) {
                    urms.pvalue <- 1 - pnorm(urms.z)
                }
            }
        }
    }

    # labels
    if(type == "rmr") {
        OUT <- list(rmr = rms, rmr.se = rms.se,
                    rmr.exactfit.z = rms.z, rmr.exactfit.pvalue = rms.pvalue,
                    urmr = urms, urmr.se = urms.se,
                    urmr.ci.lower = urms.ci.lower,
                    urmr.ci.upper = urms.ci.upper,
                    urmr.closefit.h0.value = unbiased.test.val,
                    urmr.closefit.z = urms.z,
                    urmr.closefit.pvalue = urms.pvalue)
    } else if(type == "srmr") {
        OUT <- list(srmr = rms, srmr.se = rms.se,
                    srmr.exactfit.z = rms.z, srmr.exactfit.pvalue = rms.pvalue,
                    usrmr = urms, usrmr.se = urms.se,
                    usrmr.ci.lower = urms.ci.lower,
                    usrmr.ci.upper = urms.ci.upper,
                    usrmr.closefit.h0.value = unbiased.test.val,
                    usrmr.closefit.z = urms.z,
                    usrmr.closefit.pvalue = urms.pvalue)
    } else if(type == "crmr") {
        OUT <- list(crmr = rms, crmr.se = rms.se,
                    crmr.exactfit.z = rms.z, crmr.exactfit.pvalue = rms.pvalue,
                    ucrmr = urms, ucrmr.se = urms.se,
                    ucrmr.ci.lower = urms.ci.lower,
                    ucrmr.cilupper = urms.ci.upper,
                    ucrmr.closefit.h0.value = unbiased.test.val,
                    ucrmr.closefit.z = urms.z,
                    ucrmr.closefit.pvalue = urms.pvalue)
    }

    unlist(OUT)
}

# generate summary statistics for the residuals
lav_residuals_summary_old <- function(resList = NULL,
                                  add.class = FALSE, add.labels = FALSE) {

    # per block
    nblocks <- length(resList)

    for(b in seq_len(nblocks)) {

        # create new list, including with summary statistics interleaved
        x <- vector("list", length = 0L)
        nel <- length(resList[[b]])
        NAMES <- names(resList[[b]])

        for(el in seq_len(nel)) {

            EL <- resList[[b]][[el]]
            if(!is.null(NAMES)) {
                NAME <- NAMES[el]
            }

            if(is.character(EL)) {
                new.x <- list(EL);
                if(add.labels) {
                    names(new.x) <- "type"
                }
                x <- c(x, new.x)

            } else if(is.matrix(EL) && isSymmetric(EL)) {
                tmp <- na.omit(lav_matrix_vech(EL))
                rms <- sqrt(sum(tmp*tmp)/length(tmp))
                mabs <- mean(abs(tmp))
                tmp2 <- na.omit(lav_matrix_vech(EL, diagonal = FALSE))
                rms.nodiag <- sqrt(sum(tmp2*tmp2)/length(tmp2))
                mabs.nodiag <- mean(abs(tmp2))
                cov.summary <- c(rms, rms.nodiag, mabs, mabs.nodiag)
                if(add.labels) {
                    names(cov.summary) <-
                        c("rms", "rms.nodiag", "mabs", "mabs.nodiag")
                }
                if(add.class) {
                    class(cov.summary) <- c("lavaan.vector", "numeric")
                }
                new.x <- list(EL, cov.summary)
                if(add.labels && !is.null(NAMES)) {
                    names(new.x) <- c(NAME, paste0(NAME,".summary"))
                }
                x <- c(x, new.x)

            } else {
                tmp <- na.omit(EL)
                rms <- sqrt(sum(tmp*tmp)/length(tmp))
                mabs <- mean(abs(tmp))
                mean.summary <- c(rms, mabs)
                if(add.labels) {
                    names(mean.summary) <- c("rms", "mabs")
                }
                if(add.class) {
                    class(mean.summary) <- c("lavaan.vector", "numeric")
                }
                new.x <- list(EL, mean.summary)
                if(add.labels && !is.null(NAMES)) {
                    names(new.x) <- c(NAME, paste0(NAME,".summary"))
                }
                x <- c(x, new.x)
            }

        } # nel

        # fill in block including summary statistics
        resList[[b]] <- x
    } # nblocks

    resList
}


# x is a list with sample statistics (eg output of inspect(fit, "sampstat")
# y is another (possibly the same) list with sample statistics
#
# to avoid many 'NAs', we set the scale-factor to 1
# if the to-be-scaled value is < 1e-05 (in absolute value)
lav_residuals_rescale <- function(x, diag.cov = NULL, diag.cov2 = NULL) {

    if(is.null(diag.cov2)) {
        diag.cov2 <- diag.cov
    }

    # make sure we can take the sqrt and invert
    diag.cov[!is.finite(diag.cov)] <- NA
    diag.cov[ diag.cov < .Machine$double.eps ] <- NA
    scale.cov <- tcrossprod(1/sqrt(diag.cov))

    # for the mean, we use diag.cov2
    diag.cov2[!is.finite(diag.cov2)] <- NA
    diag.cov2[ diag.cov2 < .Machine$double.eps ] <- NA
    scale.mean <- 1/sqrt(diag.cov2)

    # rescale cov
    if(!is.null(x[["cov"]])) {
        # catch (near) zero elements in x$cov
        near.zero.idx <- which(abs(x[["cov"]]) < 1e-05)
        scale.cov[near.zero.idx] <- 1
        x[["cov"]][] <- x[["cov"]] * scale.cov
    }
    if(!is.null(x[["res.cov"]])) {
        # catch (near) zero elements in x$res.cov
        near.zero.idx <- which(abs(x[["res.cov"]]) < 1e-05)
        scale.cov[near.zero.idx] <- 1
        x[["res.cov"]][] <- x[["res.cov"]] * scale.cov
    }

    # rescale int/mean
    if(!is.null(x[["res.int"]])) {
        # catch (near) zero elements in x$res.int
        near.zero.idx <- which(abs(x[["res.int"]]) < 1e-05)
        scale.mean[near.zero.idx] <- 1
        x[["res.int"]] <- x[["res.int"]] * scale.mean
    }

    if(!is.null(x[["mean"]])) {
        # catch (near) zero elements in x$mean
        near.zero.idx <- which(abs(x[["mean"]]) < 1e-05)
        scale.mean[near.zero.idx] <- 1
        x[["mean"]] <- x[["mean"]] * scale.mean
    }

    # FIXME: do something sensible for th, slopes, ...

    x
}