R/lav_object_inspect.R

# inspect a fitted psindex object

# backward compatibility -- wrapper around lavInspect
inspect.psindex <- function(object, what = "free", ...) {
    lavInspect.psindex(object              = object,
                      what                   = what,
                      add.labels             = TRUE,
                      add.class              = TRUE,
                      drop.list.single.group = TRUE)
}

# the `tech' version: no labels, full matrices, ... for further processing
lavTech.psindex <- function(object, 
                           what                   = "free",
                           add.labels             = FALSE,
                           add.class              = FALSE,
                           list.by.group          = FALSE,  
                           drop.list.single.group = FALSE) {

    lavInspect.psindex(object, what = what,
                      add.labels = add.labels, add.class = add.class,
                      list.by.group = list.by.group,
                      drop.list.single.group =  drop.list.single.group)
}

# the `user' version: with defaults for display only
lavInspect.psindex <- function(object,
                              what                   = "free",
                              add.labels             = TRUE,
                              add.class              = TRUE,
                              list.by.group          = TRUE,
                              drop.list.single.group = TRUE) {

    # object must inherit from class psindex
    stopifnot(inherits(object, "psindex"))

    # only a single argument
    if(length(what) > 1) {
        stop("`what' arguments contains multiple arguments; only one is allowed")
    }

    # be case insensitive
    what <- tolower(what)


    #### model matrices, with different contents ####
    if(what == "free") {
        lav_object_inspect_modelmatrices(object, what = "free", 
            type = "free", add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "impute" ||
              what == "imputed") { # just to ease the transition for semTools!
        object@imputed
    } else if(what == "partable" || what == "user") {
        lav_object_inspect_modelmatrices(object, what = "free", 
            type="partable", add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "se" ||
              what == "std.err" ||
              what == "standard.errors") {
        lav_object_inspect_modelmatrices(object, what = "se",
            add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "start" || what == "starting.values") {
        lav_object_inspect_modelmatrices(object, what = "start",
            add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "est"  || what == "estimates" ||
              what == "x") {
        lav_object_inspect_modelmatrices(object, what = "est",
            add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            #list.by.group = FALSE, for semTools only
            drop.list.single.group = drop.list.single.group)
    } else if(what == "dx.free") {
        lav_object_inspect_modelmatrices(object, what = "dx.free",
            add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "dx.all") {
        lav_object_inspect_modelmatrices(object, what = "dx.all",
            add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "std" || what == "std.all" || what == "standardized") {
        lav_object_inspect_modelmatrices(object, what = "std.all",
           add.labels = add.labels, add.class = add.class,
           list.by.group = list.by.group, 
           drop.list.single.group = drop.list.single.group)
    } else if(what == "std.lv") {
        lav_object_inspect_modelmatrices(object, what = "std.lv",
           add.labels = add.labels, add.class = add.class,
           list.by.group = list.by.group, 
           drop.list.single.group = drop.list.single.group)
    } else if(what == "std.nox") {
        lav_object_inspect_modelmatrices(object, what = "std.nox",
           add.labels = add.labels, add.class = add.class,
           list.by.group = list.by.group, 
           drop.list.single.group = drop.list.single.group)


    #### parameter table ####
    } else if(what == "list") {
        parTable(object)

    #### fit indices ####
    } else if(what == "fit" ||
              what == "fitmeasures" ||
              what == "fit.measures" ||
              what == "fit.indices") {
        fitMeasures(object)


    #### modification indices ####
    } else if(what == "mi" ||
              what == "modindices" ||
              what == "modification.indices") {
        modificationIndices(object)


    #### sample statistics #####
    } else if(what == "obs" ||
              what == "observed" ||
              what == "sampstat" ||
              what == "sampstats" ||
              what == "samplestats" ||
              what == "samp" ||
              what == "sample" ||
              what == "samplestatistics") {
        lav_object_inspect_sampstat(object, h1 = FALSE,
            add.labels = add.labels, add.class = add.class, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "h1" || what == "missing.h1" || what == "sampstat.h1") {
        lav_object_inspect_sampstat(object, h1 = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)

    #### wls.est - wls.obs - wls.v ####
    } else if(what == "wls.est") {
        lav_object_inspect_wls_est(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "wls.obs") {
        lav_object_inspect_wls_obs(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "wls.v") {
        lav_object_inspect_wls_v(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)



    #### data + missingness ####
    } else if(what == "data") {
        lav_object_inspect_data(object, 
            drop.list.single.group = drop.list.single.group)
    } else if(what == "case.idx") {
        lav_object_inspect_case_idx(object,
            drop.list.single.group = drop.list.single.group) 
    } else if(what == "ngroups") {
        object@Data@ngroups
    } else if(what == "group") {
        object@Data@group
    } else if(what == "cluster") {
        object@Data@cluster
    } else if(what == "ordered") {
        object@Data@ordered
    } else if(what == "group.label") {
        object@Data@group.label
    } else if(what == "nobs") {
        unlist( object@Data@nobs )
    } else if(what == "norig") {
        unlist( object@Data@norig )
    } else if(what == "ntotal") {
        sum(unlist( object@Data@nobs ))
    } else if(what == "coverage") {
        lav_object_inspect_missing_coverage(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what %in% c("patterns", "pattern")) {
        lav_object_inspect_missing_patterns(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "empty.idx") {
        lav_object_inspect_empty_idx(object,
            drop.list.single.group = drop.list.single.group)


    #### rsquare ####
    } else if(what == "rsquare" || what == "r-square" || what == "r2") {
         lav_object_inspect_rsquare(object, 
             add.labels = add.labels, add.class = add.class,
             drop.list.single.group = drop.list.single.group)


    #### model-implied sample statistics ####
    } else if(what == "implied" || what == "fitted") {
        lav_object_inspect_implied(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "resid" || what == "res" || what == "residual" ||
              what == "residuals") {
        lav_object_inspect_residuals(object, h1 = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cov.lv" || what == "veta") {
        lav_object_inspect_cov_lv(object,
            correlation.metric = FALSE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cor.lv") {
        lav_object_inspect_cov_lv(object,
            correlation.metric = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "mean.lv" || what == "eeta") {
        lav_object_inspect_mean_lv(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cov.all") {
        lav_object_inspect_cov_all(object,
            correlation.metric = FALSE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cor.all") {
        lav_object_inspect_cov_all(object,
            correlation.metric = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cov.ov" || what == "sigma" || what == "sigma.hat") {
        lav_object_inspect_cov_ov(object,
            correlation.metric = FALSE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "cor.ov") {
        lav_object_inspect_cov_ov(object,
            correlation.metric = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "mean.ov" || what == "mu" || what == "mu.hat") {
        lav_object_inspect_mean_ov(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "th" || what == "thresholds") {
        lav_object_inspect_th(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "vy") {
        lav_object_inspect_vy(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)


    #### specific model matrices? ####
    } else if(what == "theta" || what == "theta.cov") {
        lav_object_inspect_theta(object,  correlation.metric = FALSE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "theta.cor") {
        lav_object_inspect_theta(object,  correlation.metric = TRUE,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
 


    #### convergence, meanstructure, categorical ####
    } else if(what == "converged") {
        object@optim$converged
    } else if(what == "iterations" ||
              what == "iter" ||
              what == "niter") {
        object@optim$iterations
    } else if(what == "meanstructure") {
        object@Model@meanstructure
    } else if(what == "categorical") {
        object@Model@categorical
    } else if(what == "fixed.x") {
        object@Model@fixed.x
    } else if(what == "parameterization") {
        object@Model@parameterization
    } else if(what == "npar") {
        lav_object_inspect_npar(object, type = "free")
    } else if(what == "coef") {
        # this breaks simsem and semTools -- 0.6-1 
        #lav_object_inspect_coef(object, type = "free",
        #                        add.labels = add.labels, add.class = add.class)
        lav_object_inspect_modelmatrices(object, what = "est",
            type = "free", add.labels = add.labels, add.class = add.class,
            list.by.group = list.by.group,
            drop.list.single.group = drop.list.single.group)


    #### NACOV samplestats ####
    } else if(what == "gamma") {
        lav_object_inspect_sampstat_gamma(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)


    #### gradient, Hessian, information, first.order, vcov ####
    } else if(what == "gradient") {
        lav_object_inspect_gradient(object,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "hessian") {
        lav_object_inspect_hessian(object,
            add.labels = add.labels, add.class = add.class)

    } else if(what == "information") {
        lav_object_inspect_information(object, information = "default",
            augmented = FALSE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "information.expected") {
        lav_object_inspect_information(object, information = "expected",
            augmented = FALSE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "information.observed") {
        lav_object_inspect_information(object, information = "observed",
            augmented = FALSE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "information.first.order" || 
              what == "information.firstorder"  ||
              what == "first.order") {
        lav_object_inspect_information(object, information = "first.order",
            augmented = FALSE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)

    } else if(what == "augmented.information") {
        lav_object_inspect_information(object, information = "default",
            augmented = TRUE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "augmented.information.expected") {
        lav_object_inspect_information(object, information = "expected",
            augmented = TRUE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "augmented.information.observed") {
        lav_object_inspect_information(object, information = "observed",
            augmented = TRUE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "augmented.information.first.order" ||
              what == "augmented.first.order") {
        lav_object_inspect_information(object, information = "first.order",
            augmented = TRUE, inverted = FALSE,
            add.labels = add.labels, add.class = add.class)

    } else if(what == "inverted.information") {
        lav_object_inspect_information(object, information = "default",
            augmented = TRUE, inverted = TRUE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "inverted.information.expected") {
        lav_object_inspect_information(object, information = "expected",
            augmented = TRUE, inverted = TRUE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "inverted.information.observed") {
        lav_object_inspect_information(object, information = "observed",
            augmented = TRUE, inverted = TRUE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "inverted.information.first.order" || 
              what == "inverted.first.order") {
        lav_object_inspect_information(object, information = "first.order",
            augmented = TRUE, inverted = TRUE,
            add.labels = add.labels, add.class = add.class)

    } else if(what == "vcov") {
        lav_object_inspect_vcov(object,
            standardized = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.std.all" || what == "vcov.standardized" ||
              what == "vcov.std") {
        lav_object_inspect_vcov(object,
            standardized = TRUE, type = "std.all",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.std.lv") {
        lav_object_inspect_vcov(object,
            standardized = TRUE, type = "std.lv",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.std.nox") {
        lav_object_inspect_vcov(object,
            standardized = TRUE, type = "std.nox",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.def") {
        lav_object_inspect_vcov_def(object,
            standardized = FALSE,
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.def.std.all" || what == "vcov.def.standardized" ||
              what == "vcov.def.std") {
        lav_object_inspect_vcov_def(object,
            standardized = TRUE, type = "std.all",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.def.std.lv") {
        lav_object_inspect_vcov_def(object,
            standardized = TRUE, type = "std.lv",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "vcov.def.std.nox") {
        lav_object_inspect_vcov_def(object,
            standardized = TRUE, type = "std.nox",
            add.labels = add.labels, add.class = add.class)
    } else if(what == "ugamma" || what == "ug" || what == "u.gamma") {
        lav_object_inspect_UGamma(object,
            add.labels = add.labels, add.class = add.class)

    ### jacobians ####
    } else if(what == "delta") {
        lav_object_inspect_delta(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)

    # multilevel #
    } else if(what == "icc") {
        lav_object_inspect_icc(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)
    } else if(what == "ranef") {
        lav_object_inspect_ranef(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)

    # post-checking
    } else if(what == "post.check" || what == "post") {
        lav_object_post_check(object)

    # options
    } else if(what == "options" || what == "lavoptions") {
        object@Options

    # version
    } else if(what == "version") {
        object@version

    # call
    } else if(what == "call") {
        as.list( object@call )

    # timing
    } else if(what == "timing") {
        object@timing

    # optim
    } else if(what == "optim") {
        object@optim

    # test
    } else if(what == "test") {
        object@test

    # zero cell tables
    } else if(what == "zero.cell.tables") {
        lav_object_inspect_zero_cell_tables(object,
            add.labels = add.labels, add.class = add.class,
            drop.list.single.group = drop.list.single.group)

    #### not found ####
    } else {
        stop("unknown `what' argument in inspect function: `", what, "'")
    }

}


# helper functions (mostly to deal with older 'object' that may have
# been saved somewhere)
lav_object_inspect_est <- function(object) {
    
    if(inherits(object, "psindex")) {
        # from 0.5-19, they are in the partable
        if(!is.null(object@ParTable$est)) {
            OUT <- object@ParTable$est
        } else if("Fit" %in% slotNames(object)) {
            # in < 0.5-19, we should look in @Fit@est
            OUT <- object@Fit@est
        } else {
            PT <- parTable(object)
            OUT <- rep(as.numeric(NA), length(PT$lhs))
        }
    } else {
        # try generic coef()
        OUT <- coef(object, type = "user")
        if(is.matrix(OUT)) {
            # psindexList?
            OUT <- rowMeans(OUT)
        }
    }

    OUT
}

lav_object_inspect_se <- function(object) {
    
    # from 0.5-19, they are in the partable
    if(!is.null(object@ParTable$se)) {
        OUT <- object@ParTable$se
    } else if("Fit" %in% slotNames(object)) {
        # in < 0.5-19, we should look in @Fit@se
        OUT <- object@Fit@se
    } else {
        PT <- parTable(object)
        OUT <- rep(as.numeric(NA), length(PT$lhs))
    }

    OUT
}

lav_object_inspect_start <- function(object) {

    # from 0.5-19, they are in the partable
    if(!is.null(object@ParTable$start)) {
        OUT <- object@ParTable$start
    } else {
        # in < 0.5-19, we should look in @Fit@start
        OUT <- object@Fit@start
    }

    OUT
}

lav_object_inspect_boot <- function(object) {

    # from 0.5-19. they are in a separate slot
    tmp <- try(slot(object,"boot"), silent = TRUE)
    if(inherits(tmp, "try-error")) {
        # older version of object?
        est <- lav_object_inspect_est(object)
        BOOT <- attr(est, "BOOT.COEF")
    } else {
        # 0.5-19 way
        BOOT <- object@boot$coef
    }

    BOOT
}


lav_object_inspect_modelmatrices <- function(object, what = "free",
    type = "free", add.labels = FALSE, add.class = FALSE,
    list.by.group = FALSE,
    drop.list.single.group = FALSE) {

    GLIST <- object@Model@GLIST

    if(what == "dx.free") {
        DX <- lav_model_gradient(lavmodel       = object@Model,
                                 GLIST          = NULL,
                                 lavsamplestats = object@SampleStats,
                                 lavdata        = object@Data,
                                 lavcache       = object@Cache,
                                 type           = "free",
                                 verbose        = FALSE,
                                 forcePD        = TRUE,
                                 group.weight   = TRUE,
                                 Delta          = NULL)
    } else if (what == "dx.all") {
        GLIST <- lav_model_gradient(lavmodel   = object@Model,
                                GLIST          = NULL,
                                lavsamplestats = object@SampleStats,
                                lavdata        = object@Data,
                                lavcache       = object@Cache,
                                type           = "allofthem",
                                verbose        = FALSE,
                                forcePD        = TRUE,
                                group.weight   = TRUE,
                                Delta          = NULL)
        names(GLIST) <- names(object@Model@GLIST)
    } else if (what == "std.all") {
        STD <- standardize.est.all(object)
    } else if (what == "std.lv") {
        STD <- standardize.est.lv(object)
    } else if (what == "std.nox") {
        STD <- standardize.est.all.nox(object)
    } else if(what == "se") {
        SE <- lav_object_inspect_se(object)
    } else if (what == "start") {
        START <- lav_object_inspect_start(object)
    } else if (what == "est") {
        EST <- lav_object_inspect_est(object)
    }

    for(mm in 1:length(GLIST)) {

        if(add.labels) {
            dimnames(GLIST[[mm]]) <- object@Model@dimNames[[mm]]
        }

        if(what == "free") {
            # fill in free parameter counts
            if(type == "free") {
                m.el.idx <- object@Model@m.free.idx[[mm]]
                x.el.idx <- object@Model@x.free.idx[[mm]]
            #} else if(type == "unco") {
            #    m.el.idx <- object@Model@m.unco.idx[[mm]]
            #    x.el.idx <- object@Model@x.unco.idx[[mm]]
            } else if(type == "partable") {
                m.el.idx <- object@Model@m.user.idx[[mm]]
                x.el.idx <- object@Model@x.user.idx[[mm]]
            } else {
                stop("psindex ERROR: unknown type argument:", type, )
            }
            # erase everything
            GLIST[[mm]][,] <- 0.0
            GLIST[[mm]][m.el.idx] <- x.el.idx
        } else if(what == "se") {
            # fill in standard errors
            m.user.idx <- object@Model@m.user.idx[[mm]]
            x.user.idx <- object@Model@x.user.idx[[mm]]
            # erase everything
            GLIST[[mm]][,] <- 0.0
            GLIST[[mm]][m.user.idx] <- SE[x.user.idx]
        } else if(what == "start") {
            # fill in starting values
            m.user.idx <- object@Model@m.user.idx[[mm]]
            x.user.idx <- object@Model@x.user.idx[[mm]]
            GLIST[[mm]][m.user.idx] <- START[x.user.idx]
        } else if(what == "est") {
            # fill in estimated parameter values
            m.user.idx <- object@Model@m.user.idx[[mm]]
            x.user.idx <- object@Model@x.user.idx[[mm]]
            GLIST[[mm]][m.user.idx] <- EST[x.user.idx]
        } else if(what == "dx.free") {
            # fill in derivatives free parameters
            m.el.idx <- object@Model@m.free.idx[[mm]]
            x.el.idx <- object@Model@x.free.idx[[mm]]
            # erase everything
            GLIST[[mm]][,] <- 0.0
            GLIST[[mm]][m.el.idx] <- DX[x.el.idx]
        } else if(what %in% c("std.all", "std.lv", "std.nox")) {
            m.user.idx <- object@Model@m.user.idx[[mm]]
            x.user.idx <- object@Model@x.user.idx[[mm]]
            GLIST[[mm]][m.user.idx] <- STD[x.user.idx]
        } 

        # class
        if(add.class) {
            if(object@Model@isSymmetric[mm]) {
                class(GLIST[[mm]]) <- c("psindex.matrix.symmetric", "matrix")
            } else {
                class(GLIST[[mm]]) <- c("psindex.matrix", "matrix")
            }
        }
    }

    # try to reflect `equality constraints'
    con.flag <- FALSE
    if(what == "free" && object@Model@eq.constraints) {
        # extract constraints from parameter table
        PT <- parTable(object)
        CON <-  PT[PT$op %in% c("==","<",">") ,c("lhs","op","rhs")]
        rownames(CON) <- NULL

        # replace 'labels' by parameter numbers
        ID <- lav_partable_constraints_label_id(PT)
        LABEL <- names(ID)
        for(con in 1:nrow(CON)) {
            # lhs
            LHS.labels <- all.vars(as.formula(paste("~",CON[con,"lhs"]))) 

            if(length(LHS.labels) > 0L) {
                # par id
                LHS.freeid <- ID[match(LHS.labels, LABEL)]

                # substitute
                tmp <- CON[con,"lhs"]
                for(pat in 1:length(LHS.labels)) {
                    tmp <- sub(LHS.labels[pat], LHS.freeid[pat], tmp)
                }
                CON[con,"lhs"] <- tmp
            }

            # rhs
            RHS.labels <- all.vars(as.formula(paste("~",CON[con,"rhs"])))

            if(length(RHS.labels) > 0L) {
                # par id
                RHS.freeid <- ID[match(RHS.labels, LABEL)]
                # substitute
                tmp <- CON[con,"rhs"]
                for(pat in 1:length(RHS.labels)) {
                    tmp <- sub(RHS.labels[pat], RHS.freeid[pat], tmp)
                }
                CON[con,"rhs"] <- tmp
            }
        } # con

        # add this info at the top
        #GLIST <- c(constraints = list(CON), GLIST)
        #no, not a good idea, it does not work with list.by.group
  
        # add it as a 'header' attribute?
        attr(CON, "header") <- "Note: model contains equality constraints:"
        con.flag <- TRUE
    }

    # should we group them per block?
    if(list.by.group) {
        lavsamplestats <- object@SampleStats
        lavmodel       <- object@Model
        nmat           <- lavmodel@nmat

        OUT <- vector("list", length = lavmodel@nblocks)
        for(b in seq_len(lavmodel@nblocks)) {
            # which mm belong to this block? 
            mm.in.group <- 1:nmat[b] + cumsum(c(0,nmat))[b]
            mm.names <- names( GLIST[mm.in.group] )

            OUT[[b]] <- GLIST[mm.in.group]
        }

        if(lavmodel@nblocks == 1L && drop.list.single.group) {
            OUT <- OUT[[1]]
        } else {
            if(object@Data@nlevels == 1L && 
               length(object@Data@group.label) > 0L) {
                names(OUT) <- unlist(object@Data@group.label)
            } else if(object@Data@nlevels > 1L && 
                      length(object@Data@group.label) == 0L) {
                names(OUT) <- object@Data@level.label
            }
        }
    } else {
        OUT <- GLIST
    }

    # header
    if(con.flag) {
        attr(OUT, "header") <- CON
    }

    # psindex.list
    if(add.class) {
        class(OUT) <- c("psindex.list", "list")
    }

    OUT
}




# - fixme, should we export this function?
# - since 0.5-21, conditional.x = TRUE returns residual sample statistics
#    for ML, we have both joint and residual cov/var/...; but for 
#    categorical = TRUE, we only have residual cov/var...; so, we
#    only return residual in both cases, whenever residual
# - since 0.6-1, we extract the values from the @h1 slot (if present)
lav_object_inspect_sampstat <- function(object, h1 = FALSE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    nblocks <- object@Model@nblocks
    ov.names <- object@pta$vnames$ov
    ov.names.res <- object@pta$vnames$ov.nox
    ov.names.x   <- object@pta$vnames$ov.x
    lavsamplestats <- object@SampleStats

    # if nlevels, override h1 to be TRUE
    if(object@Data@nlevels > 1L) {
        h1 <- TRUE
    }

    # check if we have a non-empty @h1 slot
    if(!"h1" %in% slotNames(object)) {
        h1 <- FALSE
    } else if(length(object@h1) == 0L) {
        h1 <- FALSE
    } else {
        H1 <- object@h1$implied
    }

    # if h1 = FALSE and nlevels > 1L, nothing can show...
    if(!h1 && object@Data@nlevels > 1L) {
        stop("psindex ERROR: sample statistics not available; refit with option h1 = TRUE")
    }

    OUT <- vector("list", length = nblocks)
    for(b in seq_len(nblocks)) {

        if(!object@Model@conditional.x) {

            # covariance matrix
            if(h1) {
                OUT[[b]]$cov  <- H1$cov[[b]]
            } else {
                OUT[[b]]$cov  <- lavsamplestats@cov[[b]]
            }
            if(add.labels && !is.null(OUT[[b]]$cov)) {
                rownames(OUT[[b]]$cov) <- colnames(OUT[[b]]$cov) <- 
                    ov.names[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$cov) <- c("psindex.matrix.symmetric", "matrix")
            }

            # mean vector
            if(h1) {
                OUT[[b]]$mean <- as.numeric(H1$mean[[b]])
            } else {
                OUT[[b]]$mean <- as.numeric(lavsamplestats@mean[[b]])
            }
            if(add.labels) {
                names(OUT[[b]]$mean) <- ov.names[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$mean) <- c("psindex.vector", "numeric")
            }

            # thresholds
            if(object@Model@categorical) {
                if(h1) {
                    OUT[[b]]$th <- as.numeric(H1$th[[b]])
                } else {
                    OUT[[b]]$th <- as.numeric(lavsamplestats@th[[b]])
                }
                if(length(object@Model@num.idx[[b]]) > 0L) {
                    NUM.idx <- which(object@Model@th.idx[[b]] == 0)
                    OUT[[b]]$th <- OUT[[b]]$th[ -NUM.idx ]
                }
                if(add.labels) {
                    names(OUT[[b]]$th) <- object@pta$vnames$th[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$th) <- c("psindex.vector", "numeric")
                }
            }
        } # !conditional.x

        else { # if conditional.x = TRUE

            # residual covariance matrix
            if(h1) {
                OUT[[b]]$res.cov  <- H1$res.cov[[b]]
            } else {
                OUT[[b]]$res.cov  <- lavsamplestats@res.cov[[b]]
            }
            if(add.labels) {
                rownames(OUT[[b]]$res.cov) <- colnames(OUT[[b]]$res.cov) <- 
                    ov.names.res[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$res.cov) <- 
                    c("psindex.matrix.symmetric", "matrix")
            }
   
            # intercepts
            if(object@Model@meanstructure) {
                if(h1) {
                    OUT[[b]]$res.int <- as.numeric(H1$res.int[[b]])
                } else {
                    OUT[[b]]$res.int <- as.numeric(lavsamplestats@res.int[[b]])
                }
                if(add.labels) {
                    names(OUT[[b]]$res.int) <- ov.names.res[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.int) <- c("psindex.vector", "numeric")
                }
            }

            # thresholds
            if(object@Model@categorical) {
                if(h1) {
                    OUT[[b]]$res.th <- as.numeric(H1$res.th[[b]])
                } else {
                    OUT[[b]]$res.th <- as.numeric(lavsamplestats@res.th[[b]])
                }
                if(length(object@Model@num.idx[[b]]) > 0L) {
                    NUM.idx <- which(object@Model@th.idx[[b]] == 0)
                    OUT[[b]]$res.th <- OUT[[b]]$res.th[ -NUM.idx ]
                }
                if(add.labels) {
                    names(OUT[[b]]$res.th) <- object@pta$vnames$th[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.th) <- c("psindex.vector", "numeric")
                }
            }

            # slopes
            if(object@Model@nexo > 0L) {
                if(h1) {
                    OUT[[b]]$res.slopes  <- H1$res.slopes[[b]]
                } else {
                    OUT[[b]]$res.slopes  <- lavsamplestats@res.slopes[[b]]
                }
                if(add.labels) {
                    rownames(OUT[[b]]$res.slopes) <- ov.names.res[[b]]
                    colnames(OUT[[b]]$res.slopes) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.slopes) <- c("psindex.matrix", "matrix")
                }
            }

            # cov.x
            if(object@Model@nexo > 0L) {
                OUT[[b]]$cov.x  <- lavsamplestats@cov.x[[b]]
                if(add.labels) {
                    rownames(OUT[[b]]$cov.x) <- ov.names.x[[b]]
                    colnames(OUT[[b]]$cov.x) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$cov.x) <- 
                        c("psindex.matrix.symmetric", "matrix")
                }
            }

            # mean.x 
            if(object@Model@nexo > 0L) {
                OUT[[b]]$mean.x <- as.numeric(object@SampleStats@mean.x[[b]])
                if(add.labels) {
                    names(OUT[[b]]$mean.x) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$mean.x) <- c("psindex.vector", "numeric")
                }
            }

        } # conditional.x

        # stochastic weights
        if(object@Model@group.w.free) {
            # to be consistent with the 'implied' values, 
            # transform so group.w is the 'log(group.freq)'
            OUT[[b]]$group.w <- 
                log(lavsamplestats@group.w[[b]] * lavsamplestats@ntotal)
            if(add.labels) {
                names(OUT[[b]]$group.w) <- "w"
            }
            if(add.class) {
                class(OUT[[b]]$group.w) <- c("psindex.vector", "numeric")
            }
        }
    }

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

    OUT
}


lav_object_inspect_data <- function(object, add.labels = FALSE,
                                    drop.list.single.group = FALSE) {

    G <- object@Data@ngroups
    OUT <- object@Data@X

    if(add.labels) {
        for(g in 1:G) {
            colnames(OUT[[g]]) <- object@Data@ov.names[[g]]
        }
    }

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT
}

lav_object_inspect_case_idx <- function(object,
                                        drop.list.single.group = FALSE) {

    G <- object@Data@ngroups
    OUT <- object@Data@case.idx

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT
}


lav_object_inspect_rsquare <- function(object, est.std.all=NULL,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    nblocks <- object@Model@nblocks
    OUT <- vector("list", length = nblocks)

    if(is.null(est.std.all)) {
        est.std.all <- standardize.est.all(object)
    }

    partable <- object@ParTable
    partable$rsquare <- 1.0 - est.std.all
    # no values > 1.0
    partable$rsquare[partable$rsquare > 1.0] <- as.numeric(NA)

    for(b in seq_len(nblocks)) {
        ind.names <- partable$rhs[ which(partable$op == "=~" & 
                                         partable$block == b) ]
        eqs.y.names <- partable$lhs[ which(partable$op == "~"  & 
                                           partable$block == b) ]
        y.names <- unique( c(ind.names, eqs.y.names) )

        idx <- which(partable$op == "~~" & partable$lhs %in% y.names &
                     partable$rhs == partable$lhs & partable$block == b)
        tmp <- partable$rsquare[idx]

        if(add.labels && length(tmp) > 0L) {
            names(tmp) <- partable$lhs[idx]
        }
        if(add.class) {
            class(tmp) <- c("psindex.vector", "numeric")
        }

        OUT[[b]] <- tmp
    }

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

    OUT 
}

# model implied sample stats
lav_object_inspect_implied <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    nblocks <- object@Model@nblocks
    ov.names <- object@pta$vnames$ov
    ov.names.res <- object@pta$vnames$ov.nox
    ov.names.x   <- object@pta$vnames$ov.x
    lavimplied <- object@implied

    OUT <- vector("list", length = nblocks)
    for(b in seq_len(nblocks)) {

        if(!object@Model@conditional.x) {

            # covariance matrix
            OUT[[b]]$cov  <- lavimplied$cov[[b]]
            if(add.labels && !is.null(OUT[[b]]$cov)) {
                rownames(OUT[[b]]$cov) <- colnames(OUT[[b]]$cov) <-
                    ov.names[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$cov) <- c("psindex.matrix.symmetric", "matrix")
            }

            # mean vector
            OUT[[b]]$mean <- as.numeric(lavimplied$mean[[b]])
            if(add.labels) {
                names(OUT[[b]]$mean) <- ov.names[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$mean) <- c("psindex.vector", "numeric")
            }

            # thresholds
            if(object@Model@categorical) {
                OUT[[b]]$th <- as.numeric(lavimplied$th[[b]])
                if(length(object@Model@num.idx[[b]]) > 0L) {
                    NUM.idx <- which(object@Model@th.idx[[b]] == 0)
                    OUT[[b]]$th <- OUT[[b]]$th[ -NUM.idx ]
                }
                if(add.labels) {
                    names(OUT[[b]]$th) <- object@pta$vnames$th[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$th) <- c("psindex.vector", "numeric")
                }
            }
        } # !conditional.x

       else { # if conditional.x = TRUE

            # residual covariance matrix
            OUT[[b]]$res.cov  <- lavimplied$res.cov[[b]]
            if(add.labels) {
                rownames(OUT[[b]]$res.cov) <- colnames(OUT[[b]]$res.cov) <-
                    ov.names.res[[b]]
            }
            if(add.class) {
                class(OUT[[b]]$res.cov) <-
                    c("psindex.matrix.symmetric", "matrix")
            }

            # intercepts
            if(object@Model@meanstructure) {
                OUT[[b]]$res.int <- as.numeric(lavimplied$res.int[[b]])
                if(add.labels) {
                    names(OUT[[b]]$res.int) <- ov.names.res[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.int) <- c("psindex.vector", "numeric")
                }
            }

            # thresholds
            if(object@Model@categorical) {
                OUT[[b]]$res.th <- as.numeric(lavimplied$res.th[[b]])
                if(length(object@Model@num.idx[[b]]) > 0L) {
                    NUM.idx <- which(object@Model@th.idx[[b]] == 0)
                    OUT[[b]]$res.th <- OUT[[b]]$res.th[ -NUM.idx ]
                }
                if(add.labels) {
                    names(OUT[[b]]$res.th) <- object@pta$vnames$th[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.th) <- c("psindex.vector", "numeric")
                }
            }

            # slopes
            if(object@Model@nexo > 0L) {
                OUT[[b]]$res.slopes  <- lavimplied$res.slopes[[b]]
                if(add.labels) {
                    rownames(OUT[[b]]$res.slopes) <- ov.names.res[[b]]
                    colnames(OUT[[b]]$res.slopes) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$res.slopes) <- c("psindex.matrix", "matrix")
                }
            }

            # cov.x
            if(object@Model@nexo > 0L) {
                OUT[[b]]$cov.x  <- object@SampleStats@cov.x[[b]]
                if(add.labels) {
                    rownames(OUT[[b]]$cov.x) <- ov.names.x[[b]]
                    colnames(OUT[[b]]$cov.x) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$cov.x) <-
                        c("psindex.matrix.symmetric", "matrix")
                }
            }

            # mean.x 
            if(object@Model@nexo > 0L) {
                OUT[[b]]$mean.x  <- as.numeric(object@SampleStats@mean.x[[b]])
                if(add.labels) {
                    names(OUT[[b]]$mean.x) <- ov.names.x[[b]]
                }
                if(add.class) {
                    class(OUT[[b]]$mean.x) <- c("psindex.vector", "numeric")
                }
            }
            

        } # conditional.x

        # stochastic weights
        if(object@Model@group.w.free) {
            OUT[[b]]$group.w <- lavimplied$group.w[[b]]
            if(add.labels) {
                names(OUT[[b]]$group.w) <- "w" # somewhat redundant
            }
            if(add.class) {
                class(OUT[[b]]$group.w) <- c("psindex.vector", "numeric")
            }
        }
    }

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

    OUT
}


# residuals: _inspect_sampstat - _inspect_implied
lav_object_inspect_residuals <- function(object, h1 = TRUE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

   # unstandardized residuals
    obsList <- lav_object_inspect_sampstat(object, h1 = h1,
                                           add.labels = add.labels,
                                           add.class  = FALSE,
                                           drop.list.single.group = FALSE)
    estList <- lav_object_inspect_implied(object,
                                          add.labels = add.labels,
                                          add.class  = FALSE,
                                          drop.list.single.group = FALSE)
    # blocks
    nblocks <- length(obsList)

    resList <- vector("list", length = nblocks)
    for(b in seq_len(nblocks)) {
        if(object@Model@conditional.x) {
            if(!is.null(estList[[b]]$res.cov)) {
                resList[[b]]$res.cov <- ( obsList[[b]]$res.cov -
                                          estList[[b]]$res.cov )
                if(add.class) {
                    class(resList[[b]]$res.cov) <- 
                        c("psindex.matrix.symmetric", "matrix")
                }
            }
            if(!is.null(estList[[b]]$res.int)) {
                 resList[[b]]$res.int <- ( obsList[[b]]$res.int -
                                           estList[[b]]$res.int )
                if(add.class) {
                    class(resList[[b]]$res.int) <-
                        c("psindex.vector", "numeric")
                }
            }
            if(!is.null(estList[[b]]$res.th)) {
                resList[[b]]$res.th  <- ( obsList[[b]]$res.th  -
                                          estList[[b]]$res.th )
                if(add.class) {
                    class(resList[[b]]$res.th) <-
                        c("psindex.vector", "numeric")
                }
            }
            if(!is.null(estList[[b]]$res.slopes)) {
                resList[[b]]$res.slopes <- ( obsList[[b]]$res.slopes -
                                             estList[[b]]$res.slopes )
                if(add.class) {
                    class(resList[[b]]$res.slopes) <-
                        c("psindex.matrix", "matrix")
                }
            }
            if(!is.null(estList[[b]]$cov.x)) {
                resList[[b]]$cov.x  <- ( obsList[[b]]$cov.x  -
                                         estList[[b]]$cov.x )
                if(add.class) {
                    class(resList[[b]]$cov.x) <- 
                        c("psindex.matrix.symmetric", "matrix")
                }
            }
            if(!is.null(estList[[b]]$mean.x)) {
                 resList[[b]]$mean.x <- ( obsList[[b]]$mean.x -
                                          estList[[b]]$mean.x )
                if(add.class) {
                    class(resList[[b]]$mean.x) <- c("psindex.vector", "numeric")
                }
            }

        # unconditional
        } else {
            if(!is.null(estList[[b]]$cov)) {
                resList[[b]]$cov <- ( obsList[[b]]$cov -
                                      estList[[b]]$cov )
                if(add.class) {
                    class(resList[[b]]$cov) <- 
                        c("psindex.matrix.symmetric", "matrix")
                }
            }
            if(!is.null(estList[[b]]$mean)) {
                 resList[[b]]$mean <- ( obsList[[b]]$mean -
                                        estList[[b]]$mean )
                 if(add.class) {
                    class(resList[[b]]$mean) <-
                        c("psindex.vector", "numeric")
                 }
            }
            if(!is.null(estList[[b]]$th)) {
                resList[[b]]$th  <- ( obsList[[b]]$th  -
                                      estList[[b]]$th )
                if(add.class) {
                    class(resList[[b]]$th) <-
                        c("psindex.vector", "numeric")
                }
            }
        }

        # free group.w
        if(!is.null(estList[[b]]$group.w)) {
            resList[[b]]$group.w <- ( obsList[[b]]$group.w  - 
                                      estList[[b]]$group.w )
            if(add.class) {
                class(resList[[b]]$group.w) <-
                    c("psindex.vector", "numeric")
            }
        }
    }

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

    OUT
}


lav_object_inspect_cov_lv <- function(object, correlation.metric = FALSE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # compute lv covar
    OUT <- computeVETA(lavmodel = object@Model, remove.dummy.lv = TRUE)

    # nblocks
    nblocks <- length(OUT)

    # cor + labels + class
    for(b in seq_len(nblocks)) {

        if(correlation.metric && nrow(OUT[[b]]) > 1L) {
            # note: cov2cor fails if matrix is empty!
            OUT[[b]] <- cov2cor(OUT[[b]])
        }

        if(add.labels) {
            colnames(OUT[[b]]) <- rownames(OUT[[b]]) <- 
                object@pta$vnames$lv[[b]]
        }

        if(add.class) {
            class(OUT[[b]]) <- c("psindex.matrix.symmetric", "matrix")
        }
    }

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

    OUT
}

lav_object_inspect_mean_lv <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # compute lv means
    OUT <- computeEETA(lavmodel       = object@Model, 
                       lavsamplestats = object@SampleStats,
                       remove.dummy.lv = TRUE)

    # nblocks
    nblocks <- length(OUT)

    # ensure numeric
    OUT <- lapply(OUT, as.numeric)

    # labels + class
    for(b in seq_len(nblocks)) {
        if(add.labels && length(OUT[[b]]) > 0L) {
            names(OUT[[b]]) <- object@pta$vnames$lv.regular[[b]]
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}

lav_object_inspect_cov_all <- function(object, correlation.metric = FALSE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # compute extended model implied covariance matrix (both ov and lv)
    OUT <- computeCOV(lavmodel = object@Model, 
                      remove.dummy.lv = TRUE)

    # nblocks
    nblocks <- length(OUT)

    # cor + labels + class
    for(b in seq_len(nblocks)) {

        if(correlation.metric && nrow(OUT[[b]]) > 1L) {
            # note: cov2cor fails if matrix is empty!
            OUT[[b]] <- cov2cor(OUT[[b]])
        }

        if(add.labels) {
            NAMES <- c(object@pta$vnames$ov.model[[b]],
                       object@pta$vnames$lv.regular[[b]])
            colnames(OUT[[b]]) <- rownames(OUT[[b]]) <- NAMES
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.matrix.symmetric", "matrix")
        }
    }

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


lav_object_inspect_cov_ov <- function(object, correlation.metric = FALSE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # get model-implied covariance matrix observed
    if(object@Model@conditional.x) {
        OUT <- object@implied$res.cov
    } else {
        OUT <- object@implied$cov
    }

    # nblocks
    nblocks <- length(OUT)

    # cor + labels + class
    for(b in seq_len(nblocks)) {
 
        if(correlation.metric && nrow(OUT[[b]]) > 1L) {
            # note: cov2cor fails if matrix is empty!
            OUT[[b]] <- cov2cor(OUT[[b]])
        }

        if(add.labels) {
            colnames(OUT[[b]]) <- rownames(OUT[[b]]) <- 
                object@pta$vnames$ov.model[[b]]
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.matrix.symmetric", "matrix")
        }
    }

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

    OUT
}

lav_object_inspect_mean_ov <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # compute lv means
    if(object@Model@conditional.x) {
        OUT <- object@implied$res.int
    } else {
        OUT <- object@implied$mean
    }

    # nblocks
    nblocks <- length(OUT)
   
    # make numeric
    OUT <- lapply(OUT, as.numeric)

    # labels + class
    for(b in seq_len(nblocks)) {
        if(add.labels && length(OUT[[b]]) > 0L) {
            names(OUT[[b]]) <- object@pta$vnames$ov.model[[b]]
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}

lav_object_inspect_th <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # thresholds
    if(object@Model@conditional.x) {
        OUT <- object@implied$res.th
    } else {
        OUT <- object@implied$th
    }

    # nblocks
    nblocks <- length(OUT)

    # make numeric
    OUT <- lapply(OUT, as.numeric)

    # labels + class
    for(b in seq_len(nblocks)) {
        if(length(object@Model@num.idx[[b]]) > 0L) {
            NUM.idx <- which(object@Model@th.idx[[b]] == 0)
            OUT[[b]] <- OUT[[b]][ -NUM.idx ]
        }
        if(add.labels && length(OUT[[b]]) > 0L) {
            names(OUT[[b]]) <- object@pta$vnames$th[[b]]
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}

lav_object_inspect_vy <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # 'unconditional' model-implied variances
    #  - same as diag(Sigma.hat) if all Y are continuous)
    #  - 1.0 (or delta^2) if categorical
    #  - if also Gamma, cov.x is used (only if categorical)

    OUT <- computeVY(lavmodel = object@Model, GLIST = NULL, 
                     diagonal.only = TRUE)

    # nblocks
    nblocks <- length(OUT)

    # labels + class
    for(b in seq_len(nblocks)) {
        if(add.labels && length(OUT[[b]]) > 0L) {
            if(object@Model@categorical) {
                 names(OUT[[b]]) <- object@pta$vnames$ov.nox[[b]]
            } else {
                 names(OUT[[b]]) <- object@pta$vnames$ov[[b]]
            }
        }
        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}


lav_object_inspect_theta <- function(object, correlation.metric = FALSE,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    # get residual covariances
    OUT <- computeTHETA(lavmodel = object@Model)

    # nblocks
    nblocks <- length(OUT)

    # labels + class
    for(b in seq_len(nblocks)) {
        
        if(correlation.metric && nrow(OUT[[b]]) > 0L) {
            OUT[[b]] <- cov2cor(OUT[[b]])
        }

        if(add.labels && length(OUT[[b]]) > 0L) {
            colnames(OUT[[b]]) <- rownames(OUT[[b]]) <- 
                object@pta$vnames$ov.model[[b]]
        }

        if(add.class) {
            class(OUT[[b]]) <- c("psindex.matrix.symmetric", "matrix")
        }
    }

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

    OUT
}


lav_object_inspect_missing_coverage <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    G <- object@Data@ngroups
    OUT <- vector("list", G)
   
    for(g in 1:G) {
        if(!is.null(object@Data@Mp[[g]])) {
            OUT[[g]] <- object@Data@Mp[[g]]$coverage
        } else {
            nvar <- length(object@Data@ov.names[[g]])
            OUT[[g]] <- matrix(1.0, nvar, nvar)
        }

        if(add.labels && length(OUT[[g]]) > 0L) {
            colnames(OUT[[g]]) <- rownames(OUT[[g]]) <-
                object@pta$vnames$ov.model[[g]]
        }

        if(add.class) {
            class(OUT[[g]]) <- c("psindex.matrix.symmetric", "matrix")
        }
    }

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT
}

lav_object_inspect_missing_patterns <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    G <- object@Data@ngroups
    OUT <- vector("list", G)
   
    for(g in 1:G) {
        if(!is.null(object@Data@Mp[[g]])) {
            OUT[[g]] <- object@Data@Mp[[g]]$pat
        } else {
            nvar <- length(object@Data@ov.names[[g]])
            OUT[[g]] <- matrix(TRUE, 1L, nvar)
            rownames(OUT[[g]]) <- object@Data@nobs[[g]]
        }

        if(add.labels && length(OUT[[g]]) > 0L) {
            colnames(OUT[[g]]) <- object@pta$vnames$ov.model[[g]]
        }

        if(add.class) {
            class(OUT[[g]]) <- c("psindex.matrix", "matrix")
        }
    }

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT
}

lav_object_inspect_empty_idx <- function(object,
                                         drop.list.single.group = FALSE) {
   
    G <- object@Data@ngroups

    # get empty idx
    OUT <- vector("list", G)

    for(g in 1:G) {
        if(!is.null(object@Data@Mp[[g]])) {
            OUT[[g]] <- object@Data@Mp[[g]]$empty.idx
        } else {
            OUT[[g]] <- integer(0L)
        }
    }
    
    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT
}


lav_object_inspect_wls_est <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    OUT <- lav_model_wls_est(object@Model)

    # nblocks
    nblocks <- length(OUT)

    for(b in seq_len(nblocks)) {
        if(add.labels && length(OUT[[b]]) > 0L) {
            #FIXME!!!!
            #names(OUT[[b]]) <- ??
        }

        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}

lav_object_inspect_wls_obs <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    OUT <- object@SampleStats@WLS.obs ### FIXME: should be in @h1??

    # nblocks
    nblocks <- length(OUT)

    for(b in seq_len(nblocks)) {
        if(add.labels && length(OUT[[b]]) > 0L) {
            #FIXME!!!!
            #names(OUT[[b]]) <- ??
        }

        if(add.class) {
            class(OUT[[b]]) <- c("psindex.vector", "numeric")
        }
    }

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

    OUT
}

lav_object_inspect_wls_v <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    #OUT <- lav_model_wls_v(lavmodel       = object@Model,
    #                       lavsamplestats = object@SampleStats,
    #                       structured     = TRUE,
    #                       lavdata        = object@Data)
    # WLS.V == (traditionally) h1 expected information
    OUT <- lav_model_h1_information_expected(lavobject = object)
    # this affects fit measures gfi, agfi, pgfi


    # nblocks
    nblocks <- length(OUT)

    # if estimator == "DWLS" or "ULS", we only stored the diagonal
    # hence, we create a full matrix here
    if(object@Options$estimator %in% c("DWLS", "ULS")) {
        OUT <- lapply(OUT, 
            function(x) { nr = NROW(x); diag(x, nrow=nr, ncol=nr) })
    }

    # label + class
    for(b in seq_len(nblocks)) {
        if(add.labels && nrow(OUT[[b]]) > 0L) {
            #FIXME!!!!
            #names(OUT[[b]]) <- ??
        }

        if(add.class) {
            class(OUT[[b]]) <- c("psindex.matrix", "matrix")
        }
    }

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

    OUT
}


lav_object_inspect_sampstat_gamma <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    if(!is.null(object@SampleStats@NACOV[[1]])) {
        OUT <- object@SampleStats@NACOV
    } else {
        OUT <- lavGamma(object)
    }

    # nblocks
    nblocks <- length(OUT)

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

    OUT
}


lav_object_inspect_gradient <- function(object,
    add.labels = FALSE, add.class = FALSE) {

    if(object@SampleStats@missing.flag ||
       object@Options$estimator == "PML") {
        group.weight <- FALSE
    } else {
        group.weight <- TRUE
    }

    OUT <- lav_model_gradient(lavmodel       = object@Model,
                              GLIST          = NULL,
                              lavsamplestats = object@SampleStats,
                              lavdata        = object@Data,
                              lavcache       = object@Cache,
                              type           = "free",
                              verbose        = FALSE,
                              group.weight   = group.weight)

    # labels
    if(add.labels) {
        names(OUT) <- lav_partable_labels(object@ParTable, type="free")
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.vector", "numeric")
    }

    OUT
}

lav_object_inspect_hessian <- function(object,
    add.labels = FALSE, add.class = FALSE) {

    OUT <- lav_model_hessian(lavmodel       = object@Model,
                             lavsamplestats = object@SampleStats,
                             lavdata        = object@Data,
                             lavcache       = object@Cache,
                             lavoptions     = object@Options,
                             group.weight   = TRUE)

    # labels
    if(add.labels) {
        colnames(OUT) <- rownames(OUT) <-
            lav_partable_labels(object@ParTable, type="free")
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

lav_object_inspect_information <- function(object, 
    information = "default", augmented = FALSE, inverted = FALSE,
    add.labels = FALSE, add.class = FALSE) {

    if(information != "default") {
        # override option
        object@Options$information <- information
    } 

    OUT <- lav_model_information(lavmodel =  object@Model,
                  lavsamplestats = object@SampleStats,
                  lavdata        = object@Data,
                  lavcache       = object@Cache,
                  lavimplied     = object@implied,
                  lavoptions     = object@Options,
                  extra          = FALSE,
                  augmented      = augmented,
                  inverted       = inverted)

    # labels
    if(add.labels) {
        NAMES <- lav_partable_labels(object@ParTable, type="free")
        if(augmented) {
            nExtra <- nrow(OUT) - length(NAMES)
            if(nExtra > 0L) {
                NAMES <- c(NAMES, paste("aug", 1:nExtra, sep=""))
            }
        }
        colnames(OUT) <- rownames(OUT) <- NAMES
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

# only to provide a direct function to the old 'getVariability()' function
lav_object_inspect_firstorder <- function(object, 
    add.labels = FALSE, add.class = FALSE) {

     B0 <- lav_model_information_firstorder(lavmodel =  object@Model,
              lavsamplestats = object@SampleStats,
              lavdata        = object@Data,
              lavcache       = object@Cache,
              lavoptions     = object@Options,
              check.pd       = FALSE,
              augmented      = FALSE,
              inverted       = FALSE)
    attr(B0, "B0.group") <- NULL
    OUT <- B0

    # labels
    if(add.labels) {
        colnames(OUT) <- rownames(OUT) <-
            lav_partable_labels(object@ParTable, type="free")
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

lav_object_inspect_vcov <- function(object, standardized = FALSE,
    type = "std.all", free.only = TRUE,
    add.labels = FALSE, add.class = FALSE, remove.duplicated = FALSE) {

    npar <- max(object@ParTable$free)
    if(object@optim$npar == 0) {
        OUT <- matrix(0,0,0)
    } else {
        # check if we already have it
        tmp <- try(slot(object, "vcov"), silent = TRUE)
        if(!inherits(tmp, "try-error") && !is.null(object@vcov$vcov)) {
            OUT <- object@vcov$vcov
        } else {
        # compute it again
            OUT <- lav_model_vcov(lavmodel       = object@Model,
                                  lavsamplestats = object@SampleStats,
                                  lavoptions     = object@Options,
                                  lavdata        = object@Data,
                                  lavcache       = object@Cache,
                                  lavimplied     = object@implied,
                                  lavh1          = object@h1
                                 )
        }
    }
   
    # strip attributes
    attr(OUT, "E.inv") <- NULL
    attr(OUT, "B0") <- NULL
    attr(OUT, "B0.group") <- NULL
    attr(OUT, "Delta") <- NULL
    attr(OUT, "WLS.V") <- NULL
    attr(OUT, "BOOT.COEF") <- NULL
    attr(OUT, "BOOT.TEST") <- NULL

    # standardized?
    if(standardized) {
        if(type == "std.lv") {
            JAC <- try(lav_func_jacobian_complex(func = standardize.est.lv.x,
                x = object@optim$x, lavobject = object), silent = TRUE)
            if(inherits(JAC, "try-error")) { # eg. pnorm()
                JAC <- lav_func_jacobian_simple(func = standardize.est.lv.x,
                    x = object@optim$x, lavobject = object)
            }
        } else if(type == "std.all") {
            JAC <- try(lav_func_jacobian_complex(func = standardize.est.all.x,
                x = object@optim$x, object = object), silent = TRUE)
            if(inherits(JAC, "try-error")) { # eg. pnorm()
                JAC <- lav_func_jacobian_simple(func = standardize.est.all.x,
                    x = object@optim$x, lavobject = object)
            }
        } else if(type == "std.nox") {
            JAC <- 
                try(lav_func_jacobian_complex(func = standardize.est.all.nox.x,
                    x = object@optim$x, lavobject = object), silent = TRUE)
            if(inherits(JAC, "try-error")) { # eg. pnorm()
                JAC <- 
                    lav_func_jacobian_simple(func = standardize.est.all.nox.x,
                        x = object@optim$x, lavobject = object)
            }
        }

        # JAC contains *all* parameters in the parameter table
        if(free.only) { 
            free.idx <- which(object@ParTable$free > 0L)
            JAC <- JAC[free.idx,, drop = FALSE]
        }
        OUT <- JAC %*% OUT %*% t(JAC)
    }

    # labels
    if(add.labels) {
        colnames(OUT) <- rownames(OUT) <-
            lav_partable_labels(object@ParTable, type="free")
    }

    # alias?
    if(remove.duplicated && object@Model@eq.constraints) {
        simple.flag <- lav_constraints_check_simple(object@Model)
        if(simple.flag) {
            LAB <- lav_partable_labels(object@ParTable, type="free")
            dup.flag <- duplicated(LAB)
            OUT <- OUT[!dup.flag, !dup.flag, drop = FALSE]
        } else {
            warning("psindex WARNING: alias is TRUE, but equality constraints do not appear to be simple; returning full vcov")
        }
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

lav_object_inspect_vcov_def <- function(object, standardized = FALSE,
    type = "std.all", add.labels = FALSE, add.class = FALSE) {

    lavmodel    <- object@Model
    lavpartable <- object@ParTable
    def.idx <- which(lavpartable$op == ":=")

    if(length(def.idx) == 0L) {
        return( matrix(0,0,0) )
    }

    if(standardized) {
        # compute VCOV for "free" parameters only
        VCOV <- lav_object_inspect_vcov(object, 
                                        standardized = TRUE,
                                        type = type, free.only = FALSE,
                                        add.labels = FALSE, add.class = FALSE)
        OUT <- VCOV[def.idx, def.idx, drop = FALSE]
    } else {

        # get free parameters
        x <- lav_model_get_parameters(lavmodel, type = "free")

        # bootstrap or not?
        if(!is.null(object@boot$coef)) {
            BOOT <- object@boot$coef
            BOOT.def <- apply(BOOT, 1L, lavmodel@def.function)
            if(length(def.idx) == 1L) {
                BOOT.def <- as.matrix(BOOT.def)
            } else {
                BOOT.def <- t(BOOT.def)
            }
            OUT <- cov(BOOT.def)
        } else {
            # VCOV
            VCOV <- lav_object_inspect_vcov(object,
                                            standardized = FALSE,
                                            type = type, free.only = TRUE,
                                            add.labels = FALSE, 
                                            add.class = FALSE)
    
            # regular delta method
            JAC <- try(lav_func_jacobian_complex(func = lavmodel@def.function,
                       x = x), silent=TRUE)
            if(inherits(JAC, "try-error")) { # eg. pnorm()
                JAC <- lav_func_jacobian_simple(func = lavmodel@def.function,
                                                x = x)
            }
            OUT <- JAC %*% VCOV %*% t(JAC)
        }
    }

    # labels
    if(add.labels) {
        LHS.names <- lavpartable$lhs[def.idx]
        colnames(OUT) <- rownames(OUT) <- LHS.names
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

lav_object_inspect_UGamma <- function(object,
    add.labels = FALSE, add.class = FALSE) {

    out <- lav_test_satorra_bentler(lavobject     = object,
                                    return.ugamma = TRUE)
    OUT <- out$UGamma

    # labels
    if(add.labels) {
       # colnames(OUT) <- rownames(OUT) <-
    }

    # class
    if(add.class) {
        class(OUT) <- c("psindex.matrix.symmetric", "matrix")
    }

    OUT
}

# Delta (jacobian: d samplestats / d free_parameters)
lav_object_inspect_delta <- function(object,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    lavmodel    <- object@Model
    lavdata     <- object@Data
    lavpartable <- object@ParTable
    lavpta      <- object@pta

    OUT <- lav_object_inspect_delta_internal(lavmodel = lavmodel, 
        lavdata = lavdata, lavpartable = lavpartable, lavpta = lavpta,
        add.labels = add.labels, add.class = add.class, 
        drop.list.single.group = drop.list.single.group)   

    OUT
}

lav_object_inspect_delta_internal <- function(lavmodel = NULL, lavdata  = NULL,
    lavpartable = NULL, lavpta = NULL,
    add.labels = FALSE, add.class = FALSE, drop.list.single.group = FALSE) {

    OUT <- computeDelta(lavmodel)

    categorical    <- lavmodel@categorical
    conditional.x  <- lavmodel@conditional.x
    group.w.free   <- lavmodel@group.w.free
    nvar           <- lavmodel@nvar
    num.idx        <- lavmodel@num.idx
    th.idx         <- lavmodel@th.idx
    nexo           <- lavmodel@nexo
    nblocks        <- lavmodel@nblocks

    if(add.labels) {
        PNAMES <- lav_partable_labels(lavpartable, type="free")

        if(is.null(lavpta)) {
            lavpta <- lav_partable_attributes(lavpartable)
        }

        if(lavdata@nlevels > 1L) {
            # store names per block, rbind later
            NAMES <- vector("list", length = nblocks)
        }

        for(g in 1:nblocks) {

            if(lavdata@nlevels == 1L) {
                colnames(OUT[[g]]) <- PNAMES
            }

            if(conditional.x) {
                ov.names <- lavpta$vnames$ov.nox[[g]]
            } else {
                ov.names <- lavpta$vnames$ov[[g]]
            }
            ov.names.x <- lavpta$vnames$ov.x[[g]]
            nvar <- length(ov.names)


            names.cov <- names.cor <- names.var <- character(0L)
            names.mu <- names.pi <- names.th <- character(0L)
            names.gw <- character(0L)
 
            # Sigma
            # - if continuous: vech(Sigma)
            # - if categorical: first numeric variances, then
            #tmp <- apply(expand.grid(ov.names, ov.names), 1L, 
            #             paste, collapse = "~~")

            #if(categorical) {
            #    names.cor <- tmp[lav_matrix_vech_idx(nvar, diagonal = FALSE)]
            #    names.var <- tmp[lav_matrix_diag_idx(nvar)[num.idx[[g]]]]
            #} else {
            #    names.cov <- tmp[lav_matrix_vech_idx(nvar, diagonal = TRUE)]
            #}

            # NOTE: in 0.6-1, we use the same order, but 'label' in row-wise
            # format (eg x1 ~~ x2 instead of x2 ~~ x1)
            tmp <- matrix(apply(expand.grid(ov.names, ov.names), 1L, 
                          paste, collapse = "~~"), nrow = nvar)
            if(categorical) {
                names.cor <- lav_matrix_vechru(tmp, diagonal = FALSE)
                names.var <- diag(tmp)[num.idx[[g]]]
            } else {
                names.cov <- lav_matrix_vechru(tmp, diagonal = TRUE)
            }


            # Mu
            if(!categorical && lavmodel@meanstructure) {
                names.mu <- paste(ov.names, "~1", sep = "")
            }

            # Pi
            if(conditional.x && lavmodel@nexo[g] > 0L) {
               names.pi <- apply(expand.grid(ov.names, ov.names.x), 1L, 
                                 paste, collapse = "~")
            }

            # th
            if(categorical) {
                names.th <- lavpta$vnames$th[[g]]
                # interweave numeric intercepts, if any
                if(length(num.idx[[g]]) > 0L) {
                    tmp <- character( length(th.idx[[g]]) )
                    tmp[ th.idx[[g]] > 0 ] <- names.th
                    tmp[ th.idx[[g]] == 0 ] <- paste(ov.names[ num.idx[[g]] ],
                                                     "~1", sep = "")
                    names.th <- tmp
                }
            }

            # gw
            if(group.w.free) {
                names.gw <- "w"
            }

            if(lavdata@nlevels == 1L) {
                rownames(OUT[[g]]) <- c(names.gw,
                                        names.th, names.mu, 
                                        names.pi, 
                                        names.cov, names.var, names.cor)
                # class
                if(add.class) {
                    class(OUT[[g]]) <- c("psindex.matrix", "matrix")
                }
            } else {
                NAMES[[g]] <- c(names.gw,
                                names.th, names.mu,
                                names.pi,
                                names.cov, names.var, names.cor)
            }

        } # g
    } # labels

    # multilevel
    if(lavmodel@multilevel) {
        for(g in 1:lavmodel@ngroups) {
            if(add.labels) {
                colnames(OUT[[g]]) <- PNAMES
                rownames(OUT[[g]]) <- c(NAMES[[(g-1)*2 + 1]],
                                        NAMES[[(g-1)*2 + 2]] )
            }
            if(add.class) {
                class(OUT[[g]]) <- c("psindex.matrix", "matrix")
            }
        }
    }

    if(lavmodel@ngroups == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(lavdata@group.label) > 0L) {
            names(OUT) <- unlist(lavdata@group.label)
        }
    }

    OUT
}

lav_object_inspect_zero_cell_tables <- function(object,
            add.labels = FALSE, add.class = FALSE,
            drop.list.single.group = FALSE) {

    # categorical?
    if(!object@Model@categorical) {
        warning("psindex WARNING: no categorical variables in fitted model")
        return(invisible(list()))
    }

    lavdata <- object@Data

    # create 2-way tables
    TABLE <- lavTables(object, dimension = 2L, output = "data.frame",
                       statistic = NULL)

    # select tables with empty cells
    empty.id <- TABLE$id[which(TABLE$obs.freq == 0)]
   
    
    if(length(empty.id) == 0L) {
        # only when lavInspect() is used, give message
        if(add.class) {
            cat("(There are no tables with empty cells for this fitted model)\n")
        }
        return(invisible(list()))
    } else {
        OUT <- lav_tables_cells_format(TABLE[TABLE$id %in% empty.id,],
                   lavdata = lavdata,
                   drop.list.single.group = drop.list.single.group)
    }

    OUT
}

lav_object_inspect_coef <- function(object, type = "free", 
                                    add.labels = FALSE, add.class = FALSE) {

    if(type == "user" || type == "all") {
        type <- "user"
        idx <- 1:length( object@ParTable$lhs )
    } else if(type == "free") {
        idx <- which(object@ParTable$free > 0L & !duplicated(object@ParTable$free))
    } else {
        stop("psindex ERROR: argument `type' must be one of free or user")
    }
    EST <- lav_object_inspect_est(object)
    cof <- EST[idx]

    # labels?
    if(add.labels) {
        names(cof) <- lav_partable_labels(object@ParTable, type = type)
    }

    # class
    if(add.class) {
        class(cof) <- c("psindex.vector", "numeric")
    }

    cof
}

lav_object_inspect_npar <- function(object, type = "free") {

    if(type == "free") {
        npar <- sum(object@ParTable$free > 0L & 
                    !duplicated(object@ParTable$free))
    } else {
        npar <- length(object@ParTable$lhs)
    }

    npar
}

lav_object_inspect_icc <- function(object, add.labels = FALSE, 
                                   add.class = FALSE,
                                   drop.list.single.group = FALSE) {

    lavdata <- object@Data
    G <- lavdata@ngroups
    OUT <- vector("list", G)

    # multilevel?
    if(lavdata@nlevels == 1L) {
        stop("psindex ERROR: intraclass correlation only available for clustered data")
    }

    if(length(object@h1) == 0L) {
        stop("psindex ERROR: h1 slot is of available; refit with h1 = TRUE")
    }

    # implied statistics
    implied <- object@h1$implied

    for(g in 1:G) {
        Sigma.W <- implied$cov[[  (g-1)*lavdata@nlevels + 1 ]]
        Sigma.B <- implied$cov[[  (g-1)*lavdata@nlevels + 2 ]]

        W.diag <- diag(Sigma.W)
        B.diag <- diag(Sigma.B)

        OUT[[g]] <- numeric(length(W.diag))

        ov.names.l <- lavdata@ov.names.l[[g]]
        w.idx <- which(ov.names.l[[1]] %in% ov.names.l[[2]])
        w.names <- ov.names.l[[1]][w.idx]
        b.idx <- match(w.names, ov.names.l[[2]])

        OUT[[g]][w.idx] <- B.diag[b.idx]/(W.diag[w.idx] + B.diag[b.idx])

        # label
        if(add.labels) {
            names(OUT[[g]]) <- ov.names.l[[1]]
        }

        # class
        if(add.class) {
            class(OUT[[g]]) <- c("psindex.vector", "numeric")
        }
    } # g

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT     
}

lav_object_inspect_ranef <- function(object, add.labels = FALSE, 
                                     add.class = FALSE,
                                     drop.list.single.group = FALSE) {

    lavdata <- object@Data
    lavsamplestats <- object@SampleStats

    G <- lavdata@ngroups
    OUT <- vector("list", G)

    # multilevel?
    if(lavdata@nlevels == 1L) {
        stop("psindex ERROR: random effects only available for clustered data (in the long format)")
    }

    # implied statistics
    lavimplied <- object@implied

    for(g in 1:G) {

        Lp <- lavdata@Lp[[g]]
        YLp <- lavsamplestats@YLp[[g]]

        # implied for this group
        group.idx <- (g - 1)*lavdata@nlevels + seq_len(lavdata@nlevels)
        implied.group <- lapply(lavimplied, function(x) x[group.idx])

        # random effects (=random intercepts or cluster means)
        out <- lav_mvnorm_cluster_implied22l(Lp = Lp, implied = implied.group)
        MB.j <- lav_mvnorm_cluster_em_estep_ranef(YLp = YLp, Lp = Lp,
                        sigma.w = out$sigma.w, sigma.b = out$sigma.b,
                        sigma.zz = out$sigma.zz, sigma.yz = out$sigma.yz,
                        mu.z = out$mu.z, mu.w = out$mu.w, mu.b = out$mu.b,
                        se = FALSE)
        OUT[[g]] <- MB.j

        ov.names.l <- lavdata@ov.names.l[[g]]

        # label
        if(add.labels) {
            colnames(OUT[[g]]) <- ov.names.l[[1]]
        }

        # class
        if(add.class) {
            class(OUT[[g]]) <- c("psindex.matrix", "matrix")

        }
    } # g

    if(G == 1L && drop.list.single.group) {
        OUT <- OUT[[1]]
    } else {
        if(length(object@Data@group.label) > 0L) {
            names(OUT) <- unlist(object@Data@group.label)
        }
    }

    OUT     
}
nietsnel/psindex documentation built on June 22, 2019, 10:56 p.m.