R/reliability.R

In semTools: Useful Tools for Structural Equation Modeling

### Sunthud Pornprasertmanit, Terrence D. Jorgensen, Yves Rosseel
### Last updated: 12 February 2025



## -----------------------------
## AVE()
##   average variance extracted
##   (not a reliability index)
## -----------------------------

##' Calculate average variance extracted
##'
##' Calculate average variance extracted (AVE) per factor from `lavaan` object
##'
##' The average variance extracted (AVE) can be calculated by
##'
##' \deqn{ AVE = \frac{\bold{1}^\prime
##' \textrm{diag}\left(\Lambda\Psi\Lambda^\prime\right)\bold{1}}{\bold{1}^\prime
##' \textrm{diag}\left(\hat{\Sigma}\right) \bold{1}}, }
##'
##' Note that this formula is modified from Fornell & Larcker (1981) in the case
##' that factor variances are not 1. The proposed formula from Fornell & Larcker
##' (1981) assumes that the factor variances are 1. Note that AVE will not be
##' provided for factors consisting of items with dual loadings. AVE is the
##' property of items but not the property of factors. AVE is calculated with
##' polychoric correlations when ordinal indicators are used.
##'
##' @importFrom lavaan lavInspect
##' @importFrom methods getMethod
##'
##' @param object A [lavaan::lavaan-class] or [lavaan.mi::lavaan.mi-class] object,
##'   expected to contain only exogenous common factors (i.e., a CFA model).
##'   Cross-loadings are not allowed and will result in `NA` for any factor with
##'   indicator(s) that cross-load.
##' @param obs.var `logical` indicating whether to compute AVE using
##'   observed variances in the denominator. Setting `FALSE` triggers
##'   using model-implied variances in the denominator.
##' @param omit.imps `character` vector specifying criteria for omitting
##'   imputations from pooled results.  Can include any of
##'   `c("no.conv", "no.se", "no.npd")`, the first 2 of which are the
##'   default setting, which excludes any imputations that did not
##'   converge or for which standard errors could not be computed.  The
##'   last option (`"no.npd"`) would exclude any imputations which
##'   yielded a nonpositive definite covariance matrix for observed or
##'   latent variables, which would include any "improper solutions" such
##'   as Heywood cases.  NPD solutions are not excluded by default because
##'   they are likely to occur due to sampling error, especially in small
##'   samples.  However, gross model misspecification could also cause
##'   NPD solutions, users can compare pooled results with and without
##'   this setting as a sensitivity analysis to see whether some
##'   imputations warrant further investigation.
##' @param omit.factors `character` vector naming any common factors
##'   modeled in `object` whose indicators' AVE is not of interest.
##' @param dropSingle `logical` indicating whether to exclude factors
##'   defined by a single indicator from the returned results. If `TRUE`
##'   (default), single indicators will still be included in the `total`
##'   column when `return.total = TRUE`.
##' @param return.df `logical` indicating whether to return reliability
##'   coefficients in a `data.frame` (one row per group/level), which is
##'   possible when every model block includes the same factors (after excluding
##'   those in `omit.factors` and applying `dropSingle`).
##'
##' @return `numeric` vector of average variance extracted from indicators
##'   per factor.  For models with multiple "blocks" (any combination of groups
##'   and levels), vectors may be returned as columns in a `data.frame`
##'   with additional columns indicating the group/level (see `return.df=`
##'   argument description for caveat).
##'
##' @author
##' Terrence D. Jorgensen (University of Amsterdam; \email{TJorgensen314@@gmail.com})
##'
##' @references
##' Fornell, C., & Larcker, D. F. (1981). Evaluating structural equation models
##' with unobservable variables and measurement errors. *Journal of
##' Marketing Research, 18*(1), 39--50. \doi{10.2307/3151312}
##'
##' @seealso [compRelSEM()] for composite reliability estimates
##'
##' @examples
##' data(HolzingerSwineford1939)
##' HS9 <- HolzingerSwineford1939[ , c("x7","x8","x9")]
##' HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )
##' names(HSbinary) <- c("y7","y8","y9")
##' HS <- cbind(HolzingerSwineford1939, HSbinary)
##'
##' HS.model <- ' visual  =~ x1 + x2 + x3
##'               textual =~ x4 + x5 + x6
##'               speed   =~ y7 + y8 + y9 '
##'
##' fit <- cfa(HS.model, data = HS, ordered = c("y7","y8","y9"), std.lv = TRUE)
##'
##' ## works for factors with exclusively continuous OR categorical indicators
##' AVE(fit) # uses observed (or unconstrained polychoric/polyserial) by default
##' AVE(fit, obs.var = FALSE)
##'
##'
##' ## works for multigroup models and for multilevel models (and both)
##' data(Demo.twolevel)
##' ## assign clusters to arbitrary groups
##' Demo.twolevel$g <- ifelse(Demo.twolevel$cluster %% 2L, "type1", "type2")
##' model2 <- ' group: type1
##'   level: within
##'     fac =~ y1 + L2*y2 + L3*y3
##'   level: between
##'     fac =~ y1 + L2*y2 + L3*y3
##'
##' group: type2
##'   level: within
##'     fac =~ y1 + L2*y2 + L3*y3
##'   level: between
##'     fac =~ y1 + L2*y2 + L3*y3
##' '
##' fit2 <- sem(model2, data = Demo.twolevel, cluster = "cluster", group = "g")
##' AVE(fit2)
##'
##'@export
AVE <- function(object, obs.var = TRUE, omit.imps = c("no.conv","no.se"),
                omit.factors = character(0), dropSingle = TRUE,
                return.df = TRUE) {
  ## numbers of blocks
  ngroups <- lavInspect(object, "ngroups")
  nLevels <- lavInspect(object, "nlevels")
  nblocks <- ngroups*nLevels #FIXME: always true?

  ## labels for groups
  if (ngroups > 1L) {
    group.label <- lavInspect(object, "group.label")
    blk.g.lab <- if (!length(group.label)) paste0("g", 1:ngroups) else group.label
  } else {
    group.label <- blk.g.lab <- NULL
  }
  ## labels for clusters
  if (nLevels > 1L) {
    #FIXME? lavInspect(object, "level.label") is always ==
    #       c("within", lavInspect(object, "cluster"))
    PT <- parTable(object)
    clus.label <- unique(PT$level)
    clus.label <- clus.label[which(clus.label != "")]
    clus.label <- clus.label[which(clus.label != 0)]
    blk.clus.lab <- if (is.numeric(clus.label)) {
      c("within", lavInspect(object, "cluster"))
    } else clus.label
  } else clus.label <- blk.clus.lab <- NULL
  ## labels for blocks
  if (nblocks > 1L) {
    block.label <- paste(rep(blk.g.lab, each = nLevels), blk.clus.lab,
                         sep = if (ngroups > 1L && nLevels > 1L) "_" else "")
  } else block.label <- NULL

  ## check for categorical
  anyCategorical <- lavInspect(object, "categorical")

  if (inherits(object, "lavaan")) {
    PHI   <- lavInspect(object, "cov.lv") # common-factor variance
    EST   <- lavInspect(object, "est")    # to extract loadings
    SIGMA <- lavInspect(object,           # total variance
                        ifelse(obs.var, "sampstat", "fitted"))
    if (nblocks == 1L) {
      PHI    <- list(PHI)
      LAMBDA <- list(EST$lambda)
      SIGMA  <- list(SIGMA$cov)
    } else {
      LAMBDA <- sapply(EST, "[[", i = "lambda", simplify = FALSE)
      SIGMA <- sapply(SIGMA, "[[", i = "cov", simplify = FALSE)
    }

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    ## common-factor variance
    phiList <- object@phiList[useImps]
    if (nblocks == 1L) for (i in 1:m) phiList[[i]] <- list(phiList[[i]])
    PHI <- list()
    for (b in 1:nblocks) {
      PHI[[b]] <- Reduce("+", lapply(phiList, "[[", i = b) ) / m
    }

    ## loadings
    LAMBDA <- vector("list", nblocks)
    if (nblocks == 1L) {
      lamList <- lapply(object@coefList[useImps], "[[", i = "lambda")
      LAMBDA[[1]] <- Reduce("+", lamList) / length(lamList)
    } else {
      for (b in 1:nblocks) {
        lamList <- lapply(object@coefList[useImps], function(i) i[[b]]$lambda)
        LAMBDA[[b]] <- Reduce("+", lamList) / length(lamList)
      }
    }

    ## total variance
    if (obs.var) {
      SIGMA <- vector("list", nblocks)
      ## loop over blocks to pool saturated-model (observed) matrices
      for (b in 1:nblocks) {
        covList <- lapply(object@h1List[useImps], function(i) i$implied$cov[[b]])
        SIGMA[[b]] <- Reduce("+", covList) / m
        rownames(SIGMA[[b]]) <- colnames(SIGMA[[b]]) <- lavNames(object, block = b)
      }
    } else {
      ## pooled model-implied matrices
      if (nblocks == 1L) {
        SIGMA <- getMethod("fitted", class(object))(object)["cov"] # retain list format
      } else {
        SIGMA <- sapply(getMethod("fitted", class(object))(object),
                        "[[", "cov", simplify = FALSE)
      }
    }

  } # end lavaan vs. lavaan.mi conditional

  ## scale polychoric/polyserial to modeled LRV scale
  if (anyCategorical) {
    SDs <- sapply(getScales(object, omit.imps = omit.imps),
                  FUN = function(x) diag(1 / as.numeric(x)),
                  simplify = FALSE)

    for (b in 1:nblocks) {
      dimnames(SDs[[b]]) <- dimnames(SIGMA[[b]])
      SIGMA[[b]] <- SDs[[b]] %*% SIGMA[[b]] %*% SDs[[b]]
    }
  }

  avevar <- list()
  for (b in 1:nblocks) {
    ## extract factor and indicator names
    LY <- LAMBDA[[b]]
    allIndNames <- rownames(LY)
    allFacNames <- colnames(LY)
    myFacNames <- setdiff(allFacNames, omit.factors)
    if (dropSingle) {
      multInd <- sapply(myFacNames, function(fn) sum(LY[,fn] != 0) > 1L)
      myFacNames <- myFacNames[multInd]
    }
    subLY <- LY[ , myFacNames, drop = FALSE]
    myIndNames <- rownames(subLY)[apply(subLY, 1L, function(x) any(x != 0))]

    ## check for cross-loadings
    Xload <- apply(subLY, 1L, function(x) sum(round(x, 5) != 0) > 1L)

    avevar[[b]] <- setNames(rep(NA, length(myFacNames)), nm = myFacNames)

    ## loop over factors
    for (fn in myFacNames) {
      idx <- which(subLY[,fn] != 0)
      if (any(Xload[idx])) next # cross-loading violates AVE definition
      commonVar <- sum(subLY[idx, fn]^2) * PHI[[b]][fn, fn]
      avevar[[b]][fn] <- commonVar / sum(diag(SIGMA[[b]])[ myIndNames[idx] ])
    }

  }

  ## drop list structure?
  if (nblocks == 1L) {
    avevar <- avevar[[1]]
    class(avevar) <- c("lavaan.vector","numeric")
    return(avevar)

  } else {
    facList <- lapply(avevar, names)
    sameNames <- all(sapply(2:nblocks, function(i) {
      isTRUE(all.equal(facList[[1]], facList[[i]]))
    } ))
    if (!(sameNames && return.df)) {
      ## can't simplify, return as a list
      for (i in seq_along(avevar)) class(avevar[[i]]) <- c("lavaan.vector","numeric")
      names(avevar) <- block.label
      return(avevar)
    }
  }

  ## concatenate each factor's AVE across blocks
  facRel <- sapply(facList[[1]], simplify = FALSE, FUN = function(nn) {
    sapply(avevar, "[[", i = nn, USE.NAMES = FALSE) # extract AVE for factor i
  })
  if (ngroups > 1L && nLevels > 1L) {
    out <- data.frame(group = rep(blk.g.lab, each = nLevels),
                      level = rep(blk.clus.lab, times = ngroups),
                      facRel)
  } else if (ngroups > 1L) {
    out <- data.frame(group = blk.g.lab, facRel)
  } else if (nLevels > 1L) {
    out <- data.frame(level = blk.clus.lab, facRel)
  }
  class(out) <- c("lavaan.data.frame","data.frame")

  out
}



## ------------
## compRelSEM()
## ------------

##' Composite Reliability using SEM
##'
##' Calculate composite reliability from estimated factor-model parameters
##'
##' Several coefficients for factor-analysis reliability have been termed
##' "omega", which Cho (2021) argues is a misleading misnomer and argues for
##' using \eqn{\rho} to represent them all, differentiated by descriptive
##' subscripts.  In our package, we strive to provide unlabeled coefficients,
##' leaving it to the user to decide on a label in their report.  But we do
##' use the symbols \eqn{\alpha} and \eqn{\omega} in the formulas below in order
##' to distinguish coefficients that do (not) assume essential tau-equivalence.
##' For higher-order constructs with latent indicators, only \eqn{\omega} is
##' available. Lai's (2021) multilevel coefficients are labeled in accordance
##' with the symbols used in that article (more details below).
##'
##' Bentler (1968) first introduced factor-analysis reliability for a
##' unidimensional factor model with congeneric indicators, labeling the
##' coeficients \eqn{\alpha}.  McDonald (1999) later referred to this
##' *and other reliability coefficients*, first as \eqn{\theta} (in 1970),
##' then as \eqn{\omega}, which is a source of confusion when reporting
##' coefficients (Cho, 2021).  Coefficients based on factor models were later
##' generalized to account for multidimenisionality (possibly with
##' cross-loadings) and correlated errors. The general \eqn{\omega} formula
##' implemented in this function is:
##'
##' \deqn{ \omega = \frac{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2}
##' Var\left( \psi \right)}{\bold{1}^\prime \hat{\Sigma} \bold{1}}, }
##'
##' where \eqn{\hat{\Sigma}} can be the model-implied covariance matrix from
##' either the saturated model (i.e., the "observed" covariance matrix, used by
##' default) or from the hypothesized CFA model, controlled by the
##' `obs.var` argument. A \eqn{k}-dimensional vector \eqn{\bold{1}} is used
##' to sum elements in the matrix. Note that if the model includes any directed
##' effects (latent regression slopes), all coefficients are calculated
##' from **total** factor variances: `lavInspect(object, "cov.lv")`.
##'
##' Assuming (essential) tau-equivalence (`tau.eq=TRUE`) makes \eqn{\omega}
##' equivalent to coefficient \eqn{\alpha} from classical test theory
##' (Cronbach, 1951):
##'
##' \deqn{ \alpha = \frac{k}{k - 1}\left[ 1 - \frac{\sum^{k}_{i = 1}
##' \sigma_{ii}}{\sum^{k}_{i = 1} \sigma_{ii} + 2\sum_{i < j} \sigma_{ij}}
##' \right],}
##'
##' where \eqn{k} is the number of items in a factor's composite,
##' \eqn{\sigma_{ii}} signifies item *i*'s variance, and \eqn{\sigma_{ij}}
##' signifies the covariance between items *i* and *j*. Again, the
##' `obs.var` argument controls whether \eqn{\alpha} is calculated using
##' the observed or model-implied covariance matrix.
##'
##' By setting `return.total=TRUE`, one can estimate reliability for a
##' single composite calculated using all indicators in a multidimensional
##' CFA (Bentler, 1972, 2009). Setting `return.total = -1` will return
##' **only** the total-composite reliability (not per factor).
##'
##' **Higher-Order Factors**:
##' The reliability of a composite that represents a higher-order construct
##' requires partitioning the model-implied factor covariance matrix \eqn{\Phi}
##' in order to isolate the common-factor variance associated only with the
##' higher-order factor. Using a second-order factor model, the model-implied
##' covariance matrix of observed indicators \eqn{\hat{\Sigma}} can be
##' partitioned into 3 sources:
##' \enumerate{
##'   \item the second-order common-factor (co)variance:
##'         \eqn{\Lambda \bold{B} \Phi_2 \bold{B}^{\prime} \Lambda^{\prime}}
##'   \item the residual variance of the first-order common factors (i.e., not
##'         accounted for by the second-order factor):
##'         \eqn{\Lambda \Psi_{u} \Lambda^{\prime}}
##'   \item the measurement error of observed indicators: \eqn{\Theta}
##' }
##'
##' where \eqn{\Lambda} contains first-order factor loadings, \eqn{\bold{B}}
##' contains second-order factor loadings, \eqn{\Phi_2} is the model-implied
##' covariance matrix of the second-order factor(s), and \eqn{\Psi_{u}} is the
##' covariance matrix of first-order factor disturbances. In practice, we can
##' use the full \eqn{\bold{B}} matrix and full model-implied \eqn{\Phi} matrix
##' (i.e., including all latent factors) because the zeros in \eqn{\bold{B}}
##' will cancel out unwanted components of \eqn{\Phi}. Thus, we can calculate
##' the proportion of variance of a composite score calculated from the observed
##' indicators (e.g., a total score or scale mean) that is attributable to the
##' second-order factor (i.e., coefficient \eqn{\omega}):
##'
##' \deqn{\omega=\frac{\bold{1}^{\prime} \Lambda \bold{B} \Phi \bold{B}^{\prime}
##'   \Lambda^{\prime} \bold{1} }{ \bold{1}^{\prime} \hat{\Sigma} \bold{1}}, }
##'
##' where \eqn{\bold{1}} is the *k*-dimensional vector of 1s and *k*
##' is the number of observed indicators in the composite. Note that if a
##' higher-order factor also has observed indicators, it is necessary to model
##' the observed indicators as single-indicator constructs, so that all of the
##' higher-order factor indicators are latent (with loadings in the Beta matrix,
##' not Lambda).
##'
##' **Categorical Indicators**:
##' When all indicators (per composite) are ordinal, the `ord.scale`
##' argument controls whether the coefficient is calculated on the
##' latent-response scale (`FALSE`) or on the observed ordinal scale
##' (`TRUE`, the default).  For \eqn{\omega}-type coefficients
##' (`tau.eq=FALSE`), Green and Yang's (2009, formula 21) approach is used
##' to transform factor-model results back to the ordinal response scale. When
##' `ord.scale=TRUE` and `tau.eq=TRUE`, coefficient \eqn{\alpha} is
##' calculated using the covariance matrix calculated from the integer-valued
##' numeric weights for ordinal categories, consistent with its definition
##' (Chalmers, 2018) and the `alpha` function in the `psych` package;
##' this implies `obs.var=TRUE`, so `obs.var=FALSE` will be ignored
##' When `ord.scale=FALSE`, the standard \eqn{\alpha} formula is applied to
##' the polychoric correlation matrix ("ordinal \eqn{\alpha}"; Zumbo et al., 2007),
##' estimated from the saturated or hypothesized model (see `obs.var`),
##' and \eqn{\omega} is calculated from CFA results without applying Green and
##' Yang's (2009) correction (see Zumbo & Kroc, 2019, for a rationalization).
##' No method analogous to Green and Yang (2009) has been proposed for
##' calculating reliability with a mixture of categorical and continuous
##' indicators, so an error is returned if `object` includes factors with a
##' mixture of indicator types (unless omitted using `omit.factors`). If
##' categorical indicators load on a different factor(s) than continuous
##' indicators, then reliability will still be calculated separately for those
##' factors, but `return.total` must be `FALSE` (unless
##' `omit.factors` is used to isolate factors with indicators of the same
##' type).
##'
##' **Multilevel Measurement Models**:
##' Under the default settings, `compRelSEM()` will apply the same formula
##' in each "block" (group and/or level of analysis). In the case of multilevel
##' (ML-)SEMs, this yields "reliability" for latent within- and between-level
##' components, as proposed by Geldhof et al. (2014).  Although this works fine
##' to calculate reliability per group, this is not recommended for ML-SEMs
##' because the coefficients do not correspond to actual composites that would
##' be calculated from the observed data.  Lai (2021) proposed coefficients for
##' reliability of actual composites, depending on the type of construct, which
##' requires specifying the names of constructs for which reliability is desired
##' (or multiple constructs whose indicators would compose a multidimensional
##' composite). Configural (`config=`) and/or `shared=` constructs
##' can be specified; the same construct can be specified in both arguments, so
##' that overall scale-reliability can be estimated for a shared construct by
##' including it in `config`.  Instead of organizing the output by block
##' (the default), specifying `config=` and/or `shared=` will prompt
##' organizing the list of output by `$config` and/or `$shared`.
##'
##' \itemize{
##'   \item The overall (`_2L`) scale reliability for `config`ural
##'   constructs is returned, along with the reliability of a purely
##'   individual-level composite (`_W`, calculated by cluster-mean
##'   centering).
##'   \item The reliability for a `shared` construct quantifies
##'   generalizability across both indicators and raters (i.e., subjects rating
##'   their cluster's construct).  Lüdtke et al. (2011) refer to these as
##'   measurement error and sampling error, respectively.  An interrater
##'   reliability (IRR) coefficient is also returned, quantifying
##'   generalizability across rater/sampling-error only. To obtain a
##'   scale-reliability coefficient (quantifying a shared construct's
##'   generalizability across indicator/measurement-error only), include the
##'   same factor name in `config=`.  Jak et al. (2021) recommended
##'   modeling components of the same construct at both levels, but users may
##'   also saturate the within-level model (Lai, 2021).
##' }
##'
##' Be careful about including Level-2 variables in the model, especially
##' whether it makes sense to include them in a total composite for a Level-2
##' construct.  `dropSingle=TRUE` only prevents estimating reliability for
##' a single-indicator construct, not from including such an indicator in a
##' total composite.  It is permissible for `shared=` constructs to have
##' additional indicators at Level-2 only.  If it is necessary to model other
##' Level-2 variables (e.g., to justify the missing-at-random assumption when
##' using `missing="FIML" estimation`), they should be placed in the
##' `omit.indicators=` argument to exclude them from total composites.
##'
##'
##' @importFrom lavaan lavInspect lavNames parTable
##' @importFrom methods getMethod
##'
##' @param object A [lavaan::lavaan-class] or [lavaan.mi::lavaan.mi-class] object,
##'   expected to contain only exogenous common factors (i.e., a CFA model).
##' @param obs.var `logical` indicating whether to compute reliability
##'   using observed variances in the denominator. Setting `FALSE` triggers
##'   using model-implied variances in the denominator.
##' @param tau.eq `logical` indicating whether to assume (essential)
##'   tau-equivalence, yielding a coefficient analogous to \eqn{\alpha}.
##'   Setting `FALSE` yields an \eqn{\omega}-type coefficient.
##' @param ord.scale `logical` indicating whether to apply Green and Yang's
##'   (2009, formula 21) correction, so that reliability is calculated for the
##'   actual ordinal response scale (ignored for factors with continuous
##'   indicators).  Setting `FALSE` yields coefficients that are
##'   only applicable to the continuous latent-response scale.
##' @param config `character` vector naming any configural constructs in
##'   a multilevel CFA. For these constructs (and optional total composite),
##'   Lai's (2021) coefficients \eqn{\omega^\textrm{W}} and \eqn{\omega^\textrm{2L}}
##'   are returned (or corresponding \eqn{\alpha} coefficients when
##'   `tau.eq=TRUE`), rather than Geldhof et al.'s (2014) coefficients for
##'   hypothetical composites of latent components (although the same formula
##'   is used for \eqn{\omega^\textrm{W}} in either case). Note that the same name
##'   must be used for the factor component represented at each level of the
##'   model.
##' @param shared `character` vector naming any shared constructs in
##'   a multilevel CFA. For these constructs (and optional total composite),
##'   Lai's (2021) coefficient \eqn{\omega^\textrm{B}} or \eqn{\alpha^\textrm{B}} is
##'   returned, rather than Geldhof et al.'s (2014) between-level coefficient
##'   for hypothetical composites of latent cluster means. Lai's (2021)
##'   coefficient quantifies reliability relative to error associated with both
##'   indicators (measurement error) and subjects (sampling error), like a
##'   generalizability coefficient.  Given that subjects can be considered as
##'   raters of their cluster's shared construct, an interrater reliability
##'   (IRR) coefficient is also returned, quantifying reliability relative to
##'   rater/sampling error alone.  To quantify reliability relative to
##'   indicator/measurement error alone (i.e., \eqn{\omega^\textrm{2L}}), the
##'   `shared=` construct name(s) can additionally be included in
##'   `config=` argument.
##' @param higher `character` vector naming any higher-order constructs in
##'   `object` for which composite reliability should be calculated.
##'   Ignored when `tau.eq=TRUE` because alpha is not based on a CFA model;
##'   instead, users must fit a CFA with tau-equivalence constraints.
##'   To obtain Lai's (2021) multilevel composite-reliability indices for a
##'   higher-order factor, do not use this argument; instead, specify the
##'   higher-order factor(s) using the `shared=` or `config=` argument
##'   (`compRelSEM` will automatically check whether it includes latent
##'   indicators and apply the appropriate formula).
##' @param return.total `logical` indicating whether to return a final
##'   column containing the reliability of a composite of all indicators (not
##'   listed in `omit.indicators`) of first-order factors not listed in
##'   `omit.factors`.  Ignored in 1-factor models, and should only be set
##'   `TRUE` if all factors represent scale dimensions that could be
##'   meaningfully collapsed to a single composite (scale sum or scale mean).
##'   Setting a negative value (e.g., `-1` returns **only** the
##'   total-composite reliability (excluding coefficients per factor), except
##'   when requesting Lai's (2021) coefficients for multilevel `config`ural
##'   or `shared=` constructs.
##' @param dropSingle `logical` indicating whether to exclude factors
##'   defined by a single indicator from the returned results. If `TRUE`
##'   (default), single indicators will still be included in the `total`
##'   column when `return.total = TRUE`.
##' @param omit.factors `character` vector naming any common factors
##'   modeled in `object` whose composite reliability is not of
##'   interest. For example, higher-order or method factors. Note that
##'   [reliabilityL2()] should be used to calculate composite
##'   reliability of a higher-order factor.
##' @param omit.indicators `character` vector naming any observed variables
##'   that should be omitted from the composite whose reliability is calculated.
##' @param omit.imps `character` vector specifying criteria for omitting
##'   imputations from pooled results.  Can include any of
##'   `c("no.conv", "no.se", "no.npd")`, the first 2 of which are the
##'   default setting, which excludes any imputations that did not
##'   converge or for which standard errors could not be computed.  The
##'   last option (`"no.npd"`) would exclude any imputations which
##'   yielded a nonpositive definite covariance matrix for observed or
##'   latent variables, which would include any "improper solutions" such
##'   as Heywood cases.  NPD solutions are not excluded by default because
##'   they are likely to occur due to sampling error, especially in small
##'   samples.  However, gross model misspecification could also cause
##'   NPD solutions, users can compare pooled results with and without
##'   this setting as a sensitivity analysis to see whether some
##'   imputations warrant further investigation.
##' @param return.df `logical` indicating whether to return reliability
##'   coefficients in a `data.frame` (one row per group/level), which is
##'   possible when every model block includes the same factors (after excluding
##'   those in `omit.factors` and applying `dropSingle`).
##'
##' @return A `numeric` vector of composite reliability coefficients per
##'   factor, or a `list` of vectors per "block" (group and/or level of
##'   analysis), optionally returned as a `data.frame` when possible (see
##'   `return.df=` argument description for caveat). If there are multiple
##'   factors, whose multidimensional indicators combine into a single
##'   composite, users can request `return.total=TRUE` to add a column
##'   including a reliability coefficient for the total composite, or
##'   `return.total = -1` to return **only** the total-composite
##'   reliability (ignored when `config=` or `shared=` is specified
##'   because each factor's specification must be checked across levels).
##'
##' @author
##' Terrence D. Jorgensen (University of Amsterdam; \email{TJorgensen314@@gmail.com})
##'
##'   Uses hidden functions written by Sunthud Pornprasertmanit
##'   (\email{psunthud@@gmail.com}) for the old `reliability()` function.
##'
##' @seealso
##' [maximalRelia()] for the maximal reliability of weighted composite
##'
##' @references
##' Bentler, P. M. (1968). Alpha-maximized factor analysis (alphamax): Its
##' relation to alpha and canonical factor analysis. *Psychometrika, 33*(3),
##' 335--345. \doi{10.1007/BF02289328}
##'
##' Bentler, P. M. (1972). A lower-bound method for the dimension-free
##' measurement of internal consistency. *Social Science Research, 1*(4),
##' 343--357. \doi{10.1016/0049-089X(72)90082-8}
##'
##' Bentler, P. M. (2009). Alpha, dimension-free, and model-based internal
##' consistency reliability. *Psychometrika, 74*(1), 137--143.
##' \doi{10.1007/s11336-008-9100-1}
##'
##' Chalmers, R. P. (2018). On misconceptions and the limited usefulness of
##' ordinal alpha. *Educational and Psychological Measurement, 78*(6),
##' 1056--1071. \doi{10.1177/0013164417727036}
##'
##' Cho, E. (2021) Neither Cronbach’s alpha nor McDonald’s omega: A commentary
##' on Sijtsma and Pfadt. *Psychometrika, 86*(4), 877--886.
##' \doi{10.1007/s11336-021-09801-1}
##'
##' Cronbach, L. J. (1951). Coefficient alpha and the internal structure of
##' tests. *Psychometrika, 16*(3), 297--334. \doi{10.1007/BF02310555}
##'
##' Geldhof, G. J., Preacher, K. J., & Zyphur, M. J. (2014). Reliability
##' estimation in a multilevel confirmatory factor analysis framework.
##' *Psychological Methods, 19*(1), 72--91. \doi{10.1037/a0032138}
##'
##' Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using
##' structural equation modeling: An alternative to coefficient alpha.
##' *Psychometrika, 74*(1), 155--167. \doi{10.1007/s11336-008-9099-3}
##'
##' Jak, S., Jorgensen, T. D., & Rosseel, Y. (2021). Evaluating cluster-level
##' factor models with `lavaan` and M*plus*. *Psych, 3*(2),
##' 134--152. \doi{10.3390/psych3020012}
##'
##' Lai, M. H. C. (2021). Composite reliability of multilevel data: It’s about
##' observed scores and construct meanings. *Psychological Methods, 26*(1),
##' 90--102. \doi{10.1037/met0000287}
##'
##' Lüdtke, O., Marsh, H. W., Robitzsch, A., & Trautwein, U. (2011).
##' A 2 \eqn{\times} 2 taxonomy of multilevel latent contextual models:
##' Accuracy--bias trade-offs in full and partial error correction models.
##' *Psychological Methods, 16*(4), 444--467. \doi{10.1037/a0024376}
##'
##' McDonald, R. P. (1999). *Test theory: A unified treatment*. Mahwah, NJ:
##' Erlbaum.
##'
##' Zumbo, B. D., Gadermann, A. M., & Zeisser, C. (2007). Ordinal versions of
##' coefficients alpha and theta for Likert rating scales.
##' *Journal of Modern Applied Statistical Methods, 6*(1), 21--29.
##' \doi{10.22237/jmasm/1177992180}
##'
##' Zumbo, B. D., & Kroc, E. (2019). A measurement is a choice and Stevens’
##' scales of measurement do not help make it: A response to Chalmers.
##' *Educational and Psychological Measurement, 79*(6), 1184--1197.
##' \doi{10.1177/0013164419844305}
##'
##'
##' @examples
##'
##' data(HolzingerSwineford1939)
##' HS9 <- HolzingerSwineford1939[ , c("x7","x8","x9")]
##' HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )
##' names(HSbinary) <- c("y7","y8","y9")
##' HS <- cbind(HolzingerSwineford1939, HSbinary)
##'
##' HS.model <- ' visual  =~ x1 + x2 + x3
##'               textual =~ x4 + x5 + x6
##'               speed   =~ y7 + y8 + y9 '
##'
##' fit <- cfa(HS.model, data = HS, ordered = c("y7","y8","y9"), std.lv = TRUE)
##'
##' ## works for factors with exclusively continuous OR categorical indicators
##' compRelSEM(fit)
##'
##' ## reliability for ALL indicators only available when they are
##' ## all continuous or all categorical
##' compRelSEM(fit, omit.factors = "speed", return.total = TRUE)
##'
##'
##' ## loop over visual indicators to calculate alpha if one indicator is removed
##' for (i in paste0("x", 1:3)) {
##'   cat("Drop ", i, ":\n", sep = "")
##'   print(compRelSEM(fit, omit.factors = c("textual","speed"),
##'                    omit.indicators = i, tau.eq = TRUE))
##' }
##' ## item-total correlations obtainable by adding a composite to the data
##' HS$Visual <- HS$x1 + HS$x2 + HS$x3
##' cor(HS$Visual, y = HS[paste0("x", 1:3)])
##' ## comparable to psych::alpha(HS[paste0("x", 1:3)])
##'
##' ## Reliability of a composite that represents a higher-order factor
##' mod.hi <- ' visual  =~ x1 + x2 + x3
##'             textual =~ x4 + x5 + x6
##'             speed   =~ x7 + x8 + x9
##'             general =~ visual + textual + speed '
##'
##' fit.hi <- cfa(mod.hi, data = HolzingerSwineford1939)
##' compRelSEM(fit.hi, higher = "general")
##' ## reliabilities for lower-order composites also returned
##'
##'
##' ## works for multigroup models and for multilevel models (and both)
##' data(Demo.twolevel)
##' ## assign clusters to arbitrary groups
##' Demo.twolevel$g <- ifelse(Demo.twolevel$cluster %% 2L, "type1", "type2")
##' model2 <- ' group: type1
##'   level: 1
##'     f1 =~ y1 + L2*y2 + L3*y3
##'     f2 =~ y4 + L5*y5 + L6*y6
##'   level: 2
##'     f1 =~ y1 + L2*y2 + L3*y3
##'     f2 =~ y4 + L5*y5 + L6*y6
##'
##' group: type2
##'   level: 1
##'     f1 =~ y1 + L2*y2 + L3*y3
##'     f2 =~ y4 + L5*y5 + L6*y6
##'   level: 2
##'     f1 =~ y1 + L2*y2 + L3*y3
##'     f2 =~ y4 + L5*y5 + L6*y6
##' '
##' fit2 <- sem(model2, data = Demo.twolevel, cluster = "cluster", group = "g")
##' compRelSEM(fit2) # Geldhof's indices (hypothetical, for latent components)
##'
##' ## Lai's (2021) indices for Level-1 and configural constructs
##' compRelSEM(fit2, config = c("f1","f2"))
##' ## Lai's (2021) indices for shared (Level-2) constructs
##' ## (also an interrater reliability coefficient)
##' compRelSEM(fit2, shared = c("f1","f2"))
##'
##'
##' ## Shared construct using saturated within-level model
##' mod.sat1 <- ' level: 1
##'   y1 ~~ y1 + y2 + y3 + y4 + y5 + y6
##'   y2 ~~ y2 + y3 + y4 + y5 + y6
##'   y3 ~~ y3 + y4 + y5 + y6
##'   y4 ~~ y4 + y5 + y6
##'   y5 ~~ y5 + y6
##'   y6 ~~ y6
##'
##'   level: 2
##'   f1 =~ y1 + L2*y2 + L3*y3
##'   f2 =~ y4 + L5*y5 + L6*y6
##' '
##' fit.sat1 <- sem(mod.sat1, data = Demo.twolevel, cluster = "cluster")
##' compRelSEM(fit.sat1, shared = c("f1","f2"))
##'
##'
##' ## Simultaneous shared-and-configural model (Stapleton et al, 2016, 2019),
##' ## not recommended, but possible by omitting shared or configural factor.
##' mod.both <- ' level: 1
##'     fc =~ y1 + L2*y2 + L3*y3 + L4*y4 + L5*y5 + L6*y6
##'   level: 2
##'   ## configural construct
##'     fc =~ y1 + L2*y2 + L3*y3 + L4*y4 + L5*y5 + L6*y6
##'   ## orthogonal shared construct
##'     fs =~ NA*y1 + y2 + y3 + y4 + y5 + y6
##'     fs ~~ 1*fs + 0*fc
##' '
##' fit.both <- sem(mod.both, data = Demo.twolevel, cluster = "cluster")
##' compRelSEM(fit.both, shared = "fs", config = "fc")
##'
##' @export
compRelSEM <- function(object, obs.var = TRUE, tau.eq = FALSE, ord.scale = TRUE,
                       config = character(0), shared = character(0),
                       higher = character(0),
                       return.total = FALSE, dropSingle = TRUE,
                       omit.factors = character(0),
                       omit.indicators = character(0),
                       omit.imps = c("no.conv","no.se"), return.df = TRUE) {
  ## numbers of blocks
  ngroups <- lavInspect(object, "ngroups")
  nLevels <- lavInspect(object, "nlevels")
  nblocks <- ngroups*nLevels #FIXME: always true?
  return.total <- rep(return.total, nblocks)

  ## labels for groups
  if (ngroups > 1L) {
    group.label <- lavInspect(object, "group.label")
    blk.g.lab <- if (!length(group.label)) paste0("g", 1:ngroups) else group.label
  } else {
    group.label <- blk.g.lab <- NULL
  }
  ## labels for clusters
  if (nLevels > 1L) {
    #FIXME? lavInspect(object, "level.label") is always ==
    #       c("within", lavInspect(object, "cluster"))
    PT <- parTable(object)
    clus.label <- unique(PT$level)
    clus.label <- clus.label[which(clus.label != "")]
    clus.label <- clus.label[which(clus.label != 0)]
    blk.clus.lab <- if (is.numeric(clus.label)) {
      c("within", lavInspect(object, "cluster"))
    } else clus.label
  } else clus.label <- blk.clus.lab <- NULL
  ## labels for blocks
  if (nblocks > 1L) {
    block.label <- paste(rep(blk.g.lab, each = nLevels), blk.clus.lab,
                         sep = if (ngroups > 1L && nLevels > 1L) "_" else "")
  } else block.label <- NULL

  ## check for categorical
  anyCategorical <- lavInspect(object, "categorical")
  threshold <- if (anyCategorical) getThreshold(object, omit.imps = omit.imps) else NULL
  latScales <- if (anyCategorical) getScales(object, omit.imps = omit.imps) else NULL

  if (inherits(object, "lavaan")) {
    ## common-factor variance
    PHI <- lavInspect(object, "cov.lv") # ignored if tau.eq
    if (nblocks == 1L) PHI <- list(PHI)
    names(PHI) <- block.label

    ## factor loadings
    EST   <- lavInspect(object, "est", drop.list.single.group = FALSE)
    LAMBDA <- sapply(EST, "[[", i = "lambda", simplify = FALSE)
    names(LAMBDA) <- block.label
    ## possibly higher-order loadings?
    BETA   <- sapply(EST, "[[", i = "beta",   simplify = FALSE)
    names(BETA)   <- block.label

    ## total variance
    if (anyCategorical && tau.eq && ord.scale) {
      ## calculate SIGMA from data for alpha
      #FIXME when MLSEM available for ordinal indicators
      #     (Level 2 components available?  Extend conditional?)
      rawData <- try(lavInspect(object, "data", drop.list.single.group = FALSE),
                     silent = TRUE)
      if (inherits(rawData, "try-error"))
        stop('Error in lavInspect(fit, "data"); tau.eq= and ord.scale= cannot ',
             'both be TRUE for models fitted to summary statistics of ',
             'categorial data.')
      SIGMA <- lapply(rawData, cov)
      names(SIGMA) <- block.label

    } else {
      SIGMA <- sapply(lavInspect(object, drop.list.single.group = FALSE,
                                 what = ifelse(obs.var, "sampstat", "fitted")),
                      "[[", i = "cov", simplify = FALSE)
      names(SIGMA) <- block.label
    }

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    ## common-factor variance
    phiList <- object@phiList[useImps]
    if (nblocks == 1L) for (i in 1:m) phiList[[i]] <- list(phiList[[i]])
    PHI <- vector("list", nblocks)
    names(PHI) <- block.label
    for (b in 1:nblocks) {
      PHI[[b]] <- Reduce("+", lapply(phiList, "[[", i = b) ) / m
    }

    ## loadings (including higher-order in Beta)
    if (nblocks == 1L) {
      lamList <- lapply(object@coefList[useImps], "[[", i = "lambda")
      LAMBDA <- list(Reduce("+", lamList) / length(lamList))

      betList <- lapply(object@coefList[useImps], "[[", i = "beta")
      if (length(betList)) {
        BETA <- list(Reduce("+", betList) / length(betList))
      } else BETA <- list(NULL)

    } else {
      LAMBDA <- BETA <- vector("list", nblocks)
      names(LAMBDA) <- names(BETA) <- block.label
      for (b in 1:nblocks) {
        lamList <- lapply(object@coefList[useImps], function(i) i[[b]]$lambda)
        LAMBDA[[b]] <- Reduce("+", lamList) / length(lamList)
        betList <- lapply(object@coefList[useImps], function(i) i[[b]]$beta  )
        BETA[[b]]   <- Reduce("+", betList) / length(betList)
      }
    }

    ## total variance
    if (anyCategorical && tau.eq && ord.scale) {
      ## calculate SIGMA from data for alpha
      #FIXME when MLSEM available for ordinal indicators
      #     (Level 2 components available?  Extend conditional?)
      dataList <- object@DataList[useImps]
      SIGMA <- vector("list", nblocks)
      names(SIGMA) <- group.label #FIXME when MLSEMs can have ordinal indicators
      if (nblocks == 1L) {
        VV <- lavNames(object, type = "ov")
        impCovList <- lapply(dataList, function(DD) {
          dat <- do.call(cbind, sapply(DD[VV], as.numeric, simplify = FALSE))
          cov(dat)
        })
        SIGMA[[1]] <- Reduce("+", impCovList) / length(impCovList)

      } else {
        ## multigroup models need separate data matrices per group
        G <- lavInspect(object, "group")

        for (g in seq_along(group.label)) {
          VV <- lavNames(object, type = "ov",
                         group = ifelse(length(group.label),
                                        yes = group.label[g], no = g))

          impCovList <- lapply(dataList, function(DD) {
            RR <- DD[,G] == ifelse(length(group.label),
                                   yes = group.label[g], no = g)
            dat <- do.call(cbind, sapply(DD[RR, VV], as.numeric, simplify = FALSE))
            cov(dat)
          })
          SIGMA[[g]] <- Reduce("+", impCovList) / length(impCovList)
        } # g
      }

    } else {
      ## use model-implied SIGMA from h0 or h1 model
      if (obs.var) {
        SIGMA <- vector("list", nblocks)
        names(SIGMA) <- block.label
        ## loop over blocks to pool saturated-model (observed) matrices
        for (b in 1:nblocks) {
          covList <- lapply(object@h1List[useImps], function(i) i$implied$cov[[b]])
          SIGMA[[b]] <- Reduce("+", covList) / m
          ## The slot does not contain dimnames, so add them
          rownames(SIGMA[[b]]) <- colnames(SIGMA[[b]]) <- lavNames(object, block = b)
        }
      } else {
        ## pooled model-implied matrices
        if (nblocks == 1L) {
          SIGMA <- getMethod("fitted", class(object))(object)["cov"] # retain list format
        } else {
          SIGMA <- sapply(getMethod("fitted", class(object))(object),
                          "[[", "cov", simplify = FALSE)
          names(SIGMA) <- block.label
        }
      }

    }

  } # end lavaan vs. lavaan.mi conditional

  ## scale polychoric/polyserial to modeled LRV scale
  if (anyCategorical) {
    SDs <- sapply(latScales, function(x) diag(1 / x),
                  simplify = FALSE)
    for (b in 1:nblocks) {
      dimnames(SDs[[b]]) <- dimnames(SIGMA[[b]])
      SIGMA[[b]] <- SDs[[b]] %*% SIGMA[[b]] %*% SDs[[b]]
    }
  }

  warnTotal <- warnAlpha <- warnOmega <- FALSE
  if (!length(c(config, shared))) {

    rel <- vector("list", length = nblocks)

    for (b in 1:nblocks) {

      LY <- LAMBDA[[b]]
      allIndNames <- rownames(LY)
      allFacNames <- colnames(LY)
      myFacNames <- setdiff(allFacNames, omit.factors)
      if (dropSingle) {
        multInd <- sapply(myFacNames, function(fn) sum(LY[,fn] != 0) > 1L)
        myFacNames <- myFacNames[multInd]
      }
      subLY <- LY[ , myFacNames, drop = FALSE]
      myIndNames <- rownames(subLY)[apply(subLY, 1L, function(x) any(x != 0))]
      ## remove unwanted indicators
      myIndNames <- setdiff(myIndNames, omit.indicators)
      subLY <- subLY[myIndNames, , drop = FALSE]

      ## distinguish between categorical, continuous, and latent indicators
      nameArgs <- list(object = object)
      if (nblocks > 1L) nameArgs$block <- b
      ordNames <- do.call(lavNames, c(nameArgs, list(type = "ov.ord")))
      numNames <- do.call(lavNames, c(nameArgs, list(type = "ov.num")))
      if (anyCategorical) {
        ## identify when the (sub)set of factors are all categorical
        blockCat <- all(myIndNames %in% ordNames)
        ## identify when the (sub)set of factors have mixed indicators, so no total
        mix <- any(myIndNames %in% ordNames) && any(myIndNames %in% numNames)
      } else {
        blockCat <- FALSE
        mix <- FALSE
      }

      if (mix && return.total[b]) {
        return.total[b] <- FALSE
        if (!(tau.eq && ord.scale)) warnTotal <- TRUE
      }

      ## compute reliability per factor?
      if (return.total[b] >= 0) {

        ## set result missing by default
        rel[[b]] <- setNames(rep(NA, length(myFacNames)), nm = myFacNames)

        for (fn in myFacNames) {
          ## names of indicators with nonzero loadings
          fIndNames <- myIndNames[which(subLY[,fn] != 0)]

          ## check for ANY indicators
          if (length(fIndNames) == 0L) next
          ## check for single indicators
          if (dropSingle && length(fIndNames) == 1L) next
          ## check for categorical (or mixed) indicators
          fCat <- any(fIndNames %in% ordNames)
          ## identify when this factor has mixed indicators, so no omegas
          fMix <- fCat && any(fIndNames %in% numNames)

          ## ALPHA
          totalCov  <- SIGMA[[b]][fIndNames, fIndNames, drop = FALSE]
          if (tau.eq) {
            if (fMix && !ord.scale) {
              ## can't mix observed and latent scales
              warnAlpha <- TRUE #TODO
              next
            }
            rel[[b]][fn] <- computeAlpha(totalCov)
            next
          } # else compute omega

          ## OMEGA
          if (fMix) {
            warnOmega <- TRUE # can't (yet) mix observed and latent scales
            next
          }
          Lf <- subLY[fIndNames, fn, drop = FALSE]
          commonCov <- Lf %*% PHI[[b]][fn, fn] %*% t(Lf)
          if (fCat && ord.scale) {
            ## Green & Yang (2009)
            rel[[b]][fn] <- omegaCat(truevar = commonCov, denom = totalCov,
                                     threshold = threshold[[b]][fIndNames],
                                     scales = latScales[[b]][fIndNames])
            next
          } # else, all continuous or all LRV-scale
          rel[[b]][fn] <- sum(commonCov) / sum(totalCov)
        } # end loop over factors
      } else rel[[b]] <- c(total = as.numeric(NA))

      ## compute for total composite?
      if (return.total[b] && length(myFacNames) > 1L) {

        ## ALPHA
        totalCov  <- SIGMA[[b]][myIndNames, myIndNames, drop = FALSE]
        if (tau.eq) {
          rel[[b]]["total"] <- computeAlpha(totalCov)
          next
        } # else compute omega

        ## OMEGA
        commonCov <- subLY %*% PHI[[b]][myFacNames, myFacNames] %*% t(subLY)
        if (blockCat && ord.scale) {
          ## Green & Yang (2009)
          rel[[b]]["total"] <- omegaCat(truevar = commonCov, denom = totalCov,
                                        threshold = threshold[[b]][myIndNames],
                                        scales = latScales[[b]][myIndNames])
          next
        } # else, all continuous or all LRV-scale
        rel[[b]]["total"] <- sum(commonCov) / sum(totalCov)
      }

      ## composite(s) representing higher-order factor(s)?
      for (hf in higher) {
        ## find latent indicators
        L2 <- BETA[[b]][,hf]
        latInds <- setdiff(names(L2)[L2 != 0], omit.factors)
        ## find observed indicators
        indList <- lapply(c(hf, latInds), function(i) names(LY[,i])[ LY[,i] != 0])
        myIndNames <- setdiff(unique(do.call(c, indList)), omit.indicators)

        totalCov  <- SIGMA[[b]][myIndNames, myIndNames] # no need for drop = FALSE
        L <- LY[myIndNames, c(hf, latInds)]
        B <- BETA[[b]][c(hf, latInds), c(hf, latInds)]
        Phi <- PHI[[b]][c(hf, latInds), c(hf, latInds)]
        commonCov <- L %*% B %*% Phi %*% t(B) %*% t(L)
        if (blockCat && ord.scale) {
          ## Green & Yang (2009)
          rel[[b]][hf] <- omegaCat(truevar = commonCov, denom = totalCov,
                                   threshold = threshold[[b]][myIndNames],
                                   scales = latScales[[b]][myIndNames])
          next
        } # else, all continuous or all LRV-scale
        rel[[b]][hf] <- sum(commonCov) / sum(totalCov)

      }

    }

    ## drop list structure
    if (nblocks == 1L) {
      rel <- rel[[1]]
      class(rel) <- c("lavaan.vector","numeric")
      return(rel)

    } else {
      facList <- lapply(rel, names)
      sameNames <- all(sapply(2:nblocks, function(i) {
        isTRUE(all.equal(facList[[1]], facList[[i]]))
      } ))
      if (!(sameNames && return.df)) {
        ## can't simplify, return as a list
        for (i in seq_along(rel)) class(rel[[i]]) <- c("lavaan.vector","numeric")
        names(rel) <- block.label
        return(rel)
      }
    }

    ## concatenate each factor's reliability across blocks
    facRel <- sapply(facList[[1]], simplify = FALSE, FUN = function(nn) {
      sapply(rel, "[[", i = nn, USE.NAMES = FALSE) # extract reliability for factor i
    })
    if (ngroups > 1L && nLevels > 1L) {
      out <- data.frame(group = rep(blk.g.lab, each = nLevels),
                        level = rep(blk.clus.lab, times = ngroups),
                        facRel)
    } else if (ngroups > 1L) {
      out <- data.frame(group = blk.g.lab, facRel)
    } else if (nLevels > 1L) {
      out <- data.frame(level = blk.clus.lab, facRel)
    }
    class(out) <- c("lavaan.data.frame","data.frame")

    return(out)
  }

  if (warnTotal) {
    message('Cannot return.total when model contains both continuous and ',
            'binary/ordinal observed indicators. Use the ',
            'omit.factors= argument to choose factors with only categorical ',
            'or continuous indicators, if that is a composite of interest.\n')
  }
  if (warnAlpha) {
    message('Coefficient alpha cannot be computed for factors as though a ',
            'composite would be calculated using both observed-response scales',
            ' (for continuous indicators) and latent-response scales (for ',
            'categorical indicators).  If you want to assume tau-equivalence, ',
            'either set ord.scale=FALSE or fit a model that treats ordinal ',
            'indicators as continuous.')
  }
  if (warnOmega) {
    message('Composite reliability (omega) cannot be computed for factors ',
            'with mixed categorical and continuous indicators, unless a model ',
            'is fitted by treating ordinal indicators as continuous.')
  }


  ## otherwise, only use Lai's MULTILEVEL coefficients
  if (nLevels > 1L && length(c(config, shared))) {

    ## group-level list, each containing 2 coefs per factor/total in data.frame
    rel <- vector("list", length = ngroups)

    for (g in 1:ngroups) {

      gLab <- ifelse(length(group.label), yes = group.label[g], no = g)
      nameArgs <- list(object = object, type = "lv")
      if (ngroups > 1L) nameArgs$group <- gLab
      lv.names1 <- do.call(lavNames, c(nameArgs, list(level = clus.label[1L])))
      lv.names2 <- do.call(lavNames, c(nameArgs, list(level = clus.label[2L])))
      nameArgs$type <- "ov"
      ov.names1 <- do.call(lavNames, c(nameArgs, list(level = clus.label[1L])))
      ov.names2 <- do.call(lavNames, c(nameArgs, list(level = clus.label[2L])))
      nameArgs$type <- "lv.ind"
      allLatInds <- do.call(lavNames, c(nameArgs, list(level = clus.label[1L])))

      ## erase higher= argument, use it to collect factor names for later checks
      higher <- character(0)

      PT <- parTable(object)
      PT <- PT[PT$op == "=~", ]
      if (ngroups > 1L) PT <- PT[PT$group == gLab, ]

      ## block indices for 2 levels in this group
      #FIXME: eventually possible to model partial clustering?
      idx1 <- 1 + (g-1)*2 # within
      idx2 <- 2 + (g-1)*2 # cluster

      ## configural construct(s) defined at both levels this group?
      for (fn in config) {
        if (fn %in% omit.factors) {
          ## why would they do this?
          config <- setdiff(config, omit.factors)
          next
        }
        if (fn %in% lv.names1 && fn %in% lv.names2) {
          ## same indicators for this construct at both levels?
          indNames1 <- setdiff(PT$rhs[PT$lhs == fn & PT$level == clus.label[1L]],
                               omit.indicators)
          indNames2 <- setdiff(PT$rhs[PT$lhs == fn & PT$level == clus.label[2L]],
                               omit.indicators)
          if (!all.equal(indNames1, indNames2)) {
            stop('After removing omit.indicators=, the indicators of factor ',
                 fn, ' do not match across levels',
                 ifelse(ngroups > 1L, paste(' in group', gLab), ""))
          }
          if (dropSingle && length(indNames1) == 1L) next
          ## is it a higher-order factor?
          latInds1 <- intersect(indNames1, allLatInds)
          if (length(latInds1)) {
            lowList <- lapply(latInds1, function(fn1) {
              indNames1 <- setdiff(PT$rhs[PT$lhs == fn1 & PT$level == clus.label[1L]],
                                   omit.indicators)
              indNames2 <- setdiff(PT$rhs[PT$lhs == fn1 & PT$level == clus.label[2L]],
                                   omit.indicators)
              ## check indicators also match for lower-order factors
              if (!all.equal(indNames1, indNames2)) {
                stop('After removing omit.indicators=, the indicators of factor ',
                     fn1, ' do not match across levels',
                     ifelse(ngroups > 1L, paste(' in group', gLab), ""))
              }
              ## check indicators of lower-order factors are not also latent
              lowerLatent <- intersect(indNames1, allLatInds)
              if (length(lowerLatent))
                stop('Indicators of lower-order factor ', fn1,
                     ifelse(ngroups > 1L, paste(' in group', gLab), ""),
                     ' cannot also be latent')
              indNames1
            })
            ## update indicator list to only include relevant observed variables
            indNames1 <- intersect(c(indNames1, do.call(c, lowList)), ov.names1)
            ## update list of higher-order factors
            higher <- c(higher, fn)
          }

          Sigma1 <- SIGMA[[idx1]][indNames1, indNames1, drop = FALSE]
          Sigma2 <- SIGMA[[idx2]][indNames1, indNames1, drop = FALSE]

          if (tau.eq) {
            if (fn %in% higher) {
              warning('Cannot apply alpha (tau.eq=TRUE) to higher-order factor ',
                      fn, '. Instead, fit a higher-order CFA that imposes ',
                      'tau-equivalence constraints.')
            } else {
              ## ALPHA
              rel[[g]]$config[[fn]] <- c(`alpha_W`  = computeAlpha(Sigma1),
                                         `alpha_2L` = computeAlpha(Sigma1 + Sigma2))
            }

          } else {
            ## OMEGA
            lam1 <- LAMBDA[[idx1]][indNames1, c(fn, latInds1), drop = FALSE]
            lam2 <- LAMBDA[[idx2]][indNames1, c(fn, latInds1), drop = FALSE]
            if (!isTRUE(all.equal(lam1, lam2)))
              warning('Unequal observed-indicator loadings across levels ',
                      'detected among factors (', paste(c(fn, latInds1), collapse = ","),
                      ifelse(ngroups > 1L, paste(') in group', gLab), ")"),
                      '. omega_2L for configural constructs assumes invariance.')
            phi1 <-  PHI[[idx1]][c(fn, latInds1), c(fn, latInds1), drop = FALSE]
            phi2 <-  PHI[[idx2]][c(fn, latInds1), c(fn, latInds1), drop = FALSE]

            if (length(latInds1)) {
              bet1 <- BETA[[idx1]][c(fn, latInds1), c(fn, latInds1), drop = FALSE]
              bet2 <- BETA[[idx2]][c(fn, latInds1), c(fn, latInds1), drop = FALSE]
              if (!isTRUE(all.equal(bet1, bet2)))
                warning('Unequal higher-order loadings detected across levels ',
                        'detected for factor ', fn,
                        ifelse(ngroups > 1L, paste(' in group', gLab), ""),
                        '. omega_2L for configural constructs assumes invariance.')
              commonCov1 <- lam1 %*% bet1 %*% phi1 %*% t(bet1) %*% t(lam1)
              commonCov2 <- lam2 %*% bet2 %*% phi2 %*% t(bet2) %*% t(lam2)
            } else {
              commonCov1 <- lam1 %*% phi1 %*% t(lam1)
              commonCov2 <- lam2 %*% phi2 %*% t(lam2)
            }

            rel[[g]]$config[[fn]] <- c(`omega_W`  = sum(commonCov1)              / sum(Sigma1),
                                       `omega_2L` = sum(commonCov1 + commonCov2) / sum(Sigma1 + Sigma2))
          }

        } else {
          warning('Configural factor ', fn, 'not detected at both levels of ',
                  'analysis, so removed from config= list.  Please use the ',
                  'same name for the within- and between-level component of a ',
                  'configural construct in your syntax.')
          config <- setdiff(config, fn) # rm non-configural construct
        }
      }
      ## after removing ineligible config, still multiple for total?
      if (length(config) > 1L) {
        ## only possible if NONE of config= are higher-order factors, or if
        ## the first-order factors in config= are not latent indicators of any
        ## higher-order factors in config=
        if (return.total[idx1] && !(length(higher) && any(config %in% allLatInds))) {

          indNames1 <- setdiff(PT$rhs[PT$lhs %in% config & PT$level == clus.label[1]],
                               omit.indicators)
          ## include observed indicators of any latent indicators
          latInds1 <- intersect(indNames1, allLatInds)
          indNames1 <- setdiff(PT$rhs[PT$lhs %in% c(config, latInds1) & PT$level == clus.label[1]],
                               omit.indicators)

          Sigma1 <- SIGMA[[idx1]][indNames1, indNames1, drop = FALSE]
          Sigma2 <- SIGMA[[idx2]][indNames1, indNames1, drop = FALSE]

          if (tau.eq) {
            if (any(config %in% higher)) {
              warning('Cannot apply alpha (tau.eq=TRUE) to total composite ',
                      'that includes higher-order configural factor(s):\n',
                      paste(intersect(shared, higher), collapse = ", "),
                      '\nInstead, impose tau-equiv. in a higher-order CFA')
            } else {
              ## ALPHA
              rel[[g]]$config$total <- c(`alpha_W`  = computeAlpha(Sigma1),
                                         `alpha_2L` = computeAlpha(Sigma1 + Sigma2))
            }

          } else {
            ## OMEGA
            lam1 <- LAMBDA[[idx1]][indNames1, c(config, latInds1), drop = FALSE]
            lam2 <- LAMBDA[[idx2]][indNames1, c(config, latInds1), drop = FALSE]
            phi1 <-  PHI[[idx1]][c(config, latInds1), c(config, latInds1), drop = FALSE]
            phi2 <-  PHI[[idx2]][c(config, latInds1), c(config, latInds1), drop = FALSE]
            if (length(latInds1)) {
              bet1 <- BETA[[idx1]][c(config, latInds1), c(config, latInds1), drop = FALSE]
              bet2 <- BETA[[idx2]][c(config, latInds1), c(config, latInds1), drop = FALSE]
              commonCov1 <- lam1 %*% bet1 %*% phi1 %*% t(bet1) %*% t(lam1)
              commonCov2 <- lam2 %*% bet2 %*% phi2 %*% t(bet2) %*% t(lam2)
            } else {
              commonCov1 <- lam1 %*% phi1 %*% t(lam1)
              commonCov2 <- lam2 %*% phi2 %*% t(lam2)
            }
            rel[[g]]$config$total <- c(`omega_W`  = sum(commonCov1)              / sum(Sigma1),
                                       `omega_2L` = sum(commonCov1 + commonCov2) / sum(Sigma1 + Sigma2))
          }
        }

        ## collect 2 coefs for multiple composites into a matrix
        rel[[g]]$config <- do.call(cbind, rel[[g]]$config)
      }


      ## (reciprocal of) harmonic-mean cluster size
      Ns <- mean(1 / lavInspect(object, "cluster.size",
                                drop.list.single.group = FALSE)[[g]])
      ## reset higher= argument for later checks
      higher <- character(0)

      ## shared construct(s) defined at between level in this group?
      for (fn in shared) {
        if (fn %in% omit.factors) {
          ## why would they do this?
          shared <- setdiff(shared, omit.factors)
          next
        }
        if (fn %in% lv.names2) {
          indNames2 <- setdiff(PT$rhs[PT$lhs == fn & PT$level == clus.label[2L]],
                               omit.indicators)
          ## only Level-2 single-indicator factors are relevant to drop
          if (dropSingle && length(indNames2) == 1L) next
          ## is it a higher-order factor?
          latInds2 <- intersect(indNames2, allLatInds)
          if (length(latInds2)) {
            lowList <- lapply(latInds2, function(fn2) {
              indNames2 <- setdiff(PT$rhs[PT$lhs == fn2 & PT$level == clus.label[2L]],
                                   omit.indicators)
              ## check indicators of lower-order factors are not also latent
              lowerLatent <- intersect(indNames2, allLatInds)
              if (length(lowerLatent))
                stop('Indicators of lower-order factor ', fn2,
                     ifelse(ngroups > 1L, paste(' in group', gLab), ""),
                     ' cannot also be latent')
              indNames2
            })
            ## update indicator list to only include relevant observed variables
            indNames2 <- intersect(c(indNames2, do.call(c, lowList)), ov.names2)
            ## update list of higher-order factors
            higher <- c(higher, fn)
          }
          ## capture within-level variance components of same indicators
          ## (make sure none are Level-2 only)
          indNames1 <- intersect(indNames2, ov.names1)

          Sigma1 <- SIGMA[[idx1]][indNames1, indNames1, drop = FALSE]
          Sigma2 <- SIGMA[[idx2]][indNames2, indNames2, drop = FALSE]

          if (tau.eq) {
            if (fn %in% higher) {
              warning('Cannot apply alpha (tau.eq=TRUE) to higher-order factor ',
                      fn, '. Instead, fit a higher-order CFA that imposes ',
                      'tau-equivalence constraints.')
            } else {
              nI <- length(indNames2)
              if (nI == 1L) {
                stop('Coefficient alpha is undefined for a single indicator. ',
                     'Set tau.eq=FALSE or dropSingle=TRUE')
              }
              kw <- nI / (nI-1) # weight for alpha based on number of items
              ## ALPHA
              onlyCov2 <- Sigma2
              diag(onlyCov2) <- 0

              rel[[g]]$shared[[fn]] <- c(`alpha_B` = kw*sum(onlyCov2) / sum(Sigma1*Ns, Sigma2),
                                         `IRR`     =    sum( Sigma2 ) / sum(Sigma1*Ns, Sigma2))
            }

          } else {
            ## OMEGA
            lam2 <- LAMBDA[[idx2]][indNames2, c(fn, latInds2), drop = FALSE]
            phi2 <- PHI[[idx2]][c(fn, latInds2), c(fn, latInds2), drop = FALSE]
            if (length(latInds2)) {
              bet2 <- BETA[[idx2]][c(fn, latInds2), c(fn, latInds2), drop = FALSE]
              commonCov2 <- lam2 %*% bet2 %*% phi2 %*% t(bet2) %*% t(lam2)
            } else commonCov2 <- lam2 %*% phi2 %*% t(lam2)

            rel[[g]]$shared[[fn]] <- c(`omega_B` = sum(commonCov2) / sum(Sigma1*Ns, Sigma2),
                                       `IRR`     = sum(  Sigma2  ) / sum(Sigma1*Ns, Sigma2))
          }

        } else shared <- setdiff(shared, fn) # rm non-shared construct
      }
      ## after removing ineligible shared, still multiple for total?
      if (length(shared) > 1L) {
        ## only possible if NONE of shared= are higher-order factors, or if
        ## the first-order factors in shared= are not latent indicators of any
        ## higher-order factors in shared=
        if (return.total[idx2] && !(length(higher) && any(shared %in% allLatInds))) {
          indNames2 <- setdiff(PT$rhs[PT$lhs %in% shared & PT$level == clus.label[2]],
                               omit.indicators)
          ## include observed indicators of any latent indicators
          latInds2 <- intersect(indNames2, allLatInds)
          indNames2 <- setdiff(PT$rhs[PT$lhs %in% c(shared, latInds2) & PT$level == clus.label[2]],
                               omit.indicators)
          ## capture within-level variance components of same indicators
          ## (make sure none are Level-2 only)
          indNames1 <- intersect(indNames2, ov.names1)

          Sigma1 <- SIGMA[[idx1]][indNames1, indNames1, drop = FALSE]
          Sigma2 <- SIGMA[[idx2]][indNames2, indNames2, drop = FALSE]

          if (tau.eq) {
            if (any(shared %in% higher)) {
              warning('Cannot apply alpha (tau.eq=TRUE) to total composite ',
                      'that includes higher-order shared factor(s):\n',
                      paste(intersect(shared, higher), collapse = ", "),
                      '\nInstead, impose tau-equiv. in a higher-order CFA')
            } else {
              ## ALPHA
              nI <- length(indNames2) #TODO: justify when > length(indNames1)
                                      #     (Level-1 component exists with SD=0)
              kw <- nI / (nI-1) # weight for alpha based on number of items
              onlyCov2 <- Sigma2
              diag(onlyCov2) <- 0

              rel[[g]]$shared$total <- c(`alpha_B` = kw*sum(onlyCov2) / sum(Sigma1*Ns, Sigma2),
                                         `IRR`     =    sum( Sigma2 ) / sum(Sigma1*Ns, Sigma2))
            }

          } else {
            ## OMEGA
            lam2 <- LAMBDA[[idx2]][indNames2, c(shared, latInds2), drop = FALSE]
            phi2 <- PHI[[idx2]][c(shared, latInds2), c(shared, latInds2), drop = FALSE]
            if (length(latInds1)) {
              bet2 <- BETA[[idx2]][c(shared, latInds2), c(shared, latInds2), drop = FALSE]
              commonCov2 <- lam2 %*% bet2 %*% phi2 %*% t(bet2) %*% t(lam2)
            } else commonCov2 <- lam2 %*% phi2 %*% t(lam2)
            rel[[g]]$shared$total <- c(`omega_B` = sum(commonCov2) / sum(Sigma1*Ns, Sigma2),
                                       `IRR`     = sum(  Sigma2  ) / sum(Sigma1*Ns, Sigma2))
          }
        }

        ## collect 2 coefs for multiple composites into a matrix
        rel[[g]]$shared <- do.call(cbind, rel[[g]]$shared)
      }


    } # end loop over groups

    ## drop list structure?
    if (ngroups == 1L) {
      rel <- rel[[1]]
    } else names(rel) <- group.label

  }

  rel
}



## -------------
## reliability()
## (deprecated 10 May 2022)
## -------------


##' Composite Reliability using SEM
##'
##' Calculate composite reliability from estimated factor-model parameters
##'
##' The coefficient alpha (Cronbach, 1951) can be calculated by
##'
##' \deqn{ \alpha = \frac{k}{k - 1}\left[ 1 - \frac{\sum^{k}_{i = 1}
##' \sigma_{ii}}{\sum^{k}_{i = 1} \sigma_{ii} + 2\sum_{i < j} \sigma_{ij}}
##' \right],}
##'
##' where \eqn{k} is the number of items in a factor, \eqn{\sigma_{ii}} is the
##' item *i* observed variances, \eqn{\sigma_{ij}} is the observed
##' covariance of items *i* and *j*.
##'
##' Several coefficients for factor-analysis reliability have been termed
##' "omega", which Cho (2021) argues is a misleading misnomer and argues for
##' using \eqn{\rho} to represent them all, differentiated by descriptive
##' subscripts.  In our package, we number \eqn{\omega} based on commonly
##' applied calculations.  Bentler (1968) first introduced factor-analysis
##' reliability for a unidimensional factor model with congeneric indicators.
##' However, assuming there are no cross-loadings in a multidimensional CFA,
##' this reliability coefficient can be calculated for each factor in the model.
##'
##' \deqn{ \omega_1 =\frac{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2}
##' Var\left( \psi \right)}{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2}
##' Var\left( \psi \right) + \sum^{k}_{i = 1} \theta_{ii} + 2\sum_{i < j}
##' \theta_{ij} }, }
##'
##' where \eqn{\lambda_i} is the factor loading of item *i*, \eqn{\psi} is
##' the factor variance, \eqn{\theta_{ii}} is the variance of measurement errors
##' of item *i*, and \eqn{\theta_{ij}} is the covariance of measurement
##' errors from item *i* and *j*. McDonald (1999) later referred to
##' this *and other reliability coefficients* as "omega", which is a source
##' of confusion when reporting coefficients (Cho, 2021).
##'
##' The additional coefficients generalize the first formula by accounting for
##' multidimenisionality (possibly with cross-loadings) and correlated errors.
##' By setting `return.total=TRUE`, one can estimate reliability for a
##' single composite calculated using all indicators in the multidimensional
##' CFA (Bentler, 1972, 2009).  `"omega2"` is calculated by
##'
##' \deqn{ \omega_2 = \frac{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2}
##' Var\left( \psi \right)}{\bold{1}^\prime \hat{\Sigma} \bold{1}}, }
##'
##' where \eqn{\hat{\Sigma}} is the model-implied covariance matrix, and
##' \eqn{\bold{1}} is the \eqn{k}-dimensional vector of 1. The first and the
##' second coefficients omega will have the same value per factor in models with
##' simple structure, but they differ when there are (e.g.) cross-loadings
##' or method factors. The first coefficient omega can be viewed as the
##' reliability controlling for the other factors (like \eqn{\eta^2_{partial}} in
##' ANOVA). The second coefficient omega can be viewed as the unconditional
##' reliability (like \eqn{\eta^2} in ANOVA).
##'
##' The `"omega3"` coefficient (McDonald, 1999), sometimes referred to as
##' hierarchical omega, can be calculated by
##'
##' \deqn{ \omega_3 =\frac{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2}
##' Var\left( \psi \right)}{\bold{1}^\prime \Sigma \bold{1}}, }
##'
##' where \eqn{\Sigma} is the observed covariance matrix. If the model fits the
##' data well, \eqn{\omega_3} will be similar to \eqn{\omega_2}. Note that if
##' there is a directional effect in the model, all coefficients are calculated
##' from total factor variances: `lavInspect(object, "cov.lv")`.
##'
##' In conclusion, \eqn{\omega_1}, \eqn{\omega_2}, and \eqn{\omega_3} are
##' different in the denominator. The denominator of the first formula assumes
##' that a model is congeneric factor model where measurement errors are not
##' correlated. The second formula accounts for correlated measurement errors.
##' However, these two formulas assume that the model-implied covariance matrix
##' explains item relationships perfectly. The residuals are subject to sampling
##' error. The third formula use observed covariance matrix instead of
##' model-implied covariance matrix to calculate the observed total variance.
##' This formula is the most conservative method in calculating coefficient
##' omega.
##'
##' The average variance extracted (AVE) can be calculated by
##'
##' \deqn{ AVE = \frac{\bold{1}^\prime
##' \textrm{diag}\left(\Lambda\Psi\Lambda^\prime\right)\bold{1}}{\bold{1}^\prime
##' \textrm{diag}\left(\hat{\Sigma}\right) \bold{1}}, }
##'
##' Note that this formula is modified from Fornell & Larcker (1981) in the case
##' that factor variances are not 1. The proposed formula from Fornell & Larcker
##' (1981) assumes that the factor variances are 1. Note that AVE will not be
##' provided for factors consisting of items with dual loadings. AVE is the
##' property of items but not the property of factors. AVE is calculated with
##' polychoric correlations when ordinal indicators are used.
##'
##' Coefficient alpha is by definition applied by treating indicators as numeric
##' (see Chalmers, 2018), which is consistent with the `alpha` function in
##' the `psych` package. When indicators are ordinal, `reliability`
##' additionally applies the standard alpha calculation to the polychoric
##' correlation matrix to return Zumbo et al.'s (2007) "ordinal alpha".
##'
##' Coefficient omega for categorical items is calculated using Green and Yang's
##' (2009, formula 21) approach. Three types of coefficient omega indicate
##' different methods to calculate item total variances. The original formula
##' from Green and Yang is equivalent to \eqn{\omega_3} in this function.
##' Green and Yang did not propose a method for
##' calculating reliability with a mixture of categorical and continuous
##' indicators, and we are currently unaware of an appropriate method.
##' Therefore, when `reliability` detects both categorical and continuous
##' indicators of a factor, an error is returned. If the categorical indicators
##' load on a different factor(s) than continuous indicators, then reliability
##' will still be calculated separately for those factors, but
##' `return.total` must be `FALSE` (unless `omit.factors` is used
##' to isolate factors with indicators of the same type).
##'
##'
##' @importFrom lavaan lavInspect lavNames
##' @importFrom methods getMethod
##'
##' @param object A [lavaan::lavaan-class] or [lavaan.mi::lavaan.mi-class] object,
##'   expected to contain only exogenous common factors (i.e., a CFA model).
##' @param what `character` vector naming any reliability indices to
##'   calculate. All are returned by default. When indicators are ordinal,
##'   both traditional `"alpha"` and Zumbo et al.'s (2007) so-called
##'   "ordinal alpha" (`"alpha.ord"`) are returned, though the latter is
##'   arguably of dubious value (Chalmers, 2018).
##' @param return.total `logical` indicating whether to return a final
##'   column containing the reliability of a composite of all indicators (not
##'   listed in `omit.indicators`) of factors not listed in
##'   `omit.factors`.  Ignored in 1-factor models, and should only be set
##'   `TRUE` if all factors represent scale dimensions that could be
##'   meaningfully collapsed to a single composite (scale sum or scale mean).
##' @param dropSingle `logical` indicating whether to exclude factors
##'   defined by a single indicator from the returned results. If `TRUE`
##'   (default), single indicators will still be included in the `total`
##'   column when `return.total = TRUE`.
##' @param omit.factors `character` vector naming any common factors
##'   modeled in `object` whose composite reliability is not of
##'   interest. For example, higher-order or method factors. Note that
##'   [reliabilityL2()] should be used to calculate composite
##'   reliability of a higher-order factor.
##' @param omit.indicators `character` vector naming any observed variables
##'   that should be ignored when calculating composite reliability. This can
##'   be useful, for example, to estimate reliability when an indicator is
##'   removed.
##' @param omit.imps `character` vector specifying criteria for omitting
##'   imputations from pooled results.  Can include any of
##'   `c("no.conv", "no.se", "no.npd")`, the first 2 of which are the
##'   default setting, which excludes any imputations that did not
##'   converge or for which standard errors could not be computed.  The
##'   last option (`"no.npd"`) would exclude any imputations which
##'   yielded a nonpositive definite covariance matrix for observed or
##'   latent variables, which would include any "improper solutions" such
##'   as Heywood cases.  NPD solutions are not excluded by default because
##'   they are likely to occur due to sampling error, especially in small
##'   samples.  However, gross model misspecification could also cause
##'   NPD solutions, users can compare pooled results with and without
##'   this setting as a sensitivity analysis to see whether some
##'   imputations warrant further investigation.
##'
##' @return Reliability values (coefficient alpha, coefficients omega, average
##'   variance extracted) of each factor in each group. If there are multiple
##'   factors, a `total` column can optionally be included.
##'
##' @author Sunthud Pornprasertmanit (\email{psunthud@@gmail.com})
##'
##'   Yves Rosseel (Ghent University; \email{Yves.Rosseel@@UGent.be})
##'
##'   Terrence D. Jorgensen (University of Amsterdam; \email{TJorgensen314@@gmail.com})
##'
##' @references
##' Bentler, P. M. (1972). A lower-bound method for the dimension-free
##' measurement of internal consistency. *Social Science Research, 1*(4),
##' 343--357. \doi{10.1016/0049-089X(72)90082-8}
##'
##' Bentler, P. M. (2009). Alpha, dimension-free, and model-based internal
##' consistency reliability. *Psychometrika, 74*(1), 137--143.
##' \doi{10.1007/s11336-008-9100-1}
##'
##' Chalmers, R. P. (2018). On misconceptions and the limited usefulness of
##' ordinal alpha. *Educational and Psychological Measurement, 78*(6),
##' 1056--1071. \doi{10.1177/0013164417727036}
##'
##' Cho, E. (2021) Neither Cronbach’s alpha nor McDonald’s omega: A commentary
##' on Sijtsma and Pfadt. *Psychometrika, 86*(4), 877--886.
##' \doi{10.1007/s11336-021-09801-1}
##'
##' Cronbach, L. J. (1951). Coefficient alpha and the internal structure of
##' tests. *Psychometrika, 16*(3), 297--334. \doi{10.1007/BF02310555}
##'
##' Fornell, C., & Larcker, D. F. (1981). Evaluating structural equation models
##' with unobservable variables and measurement errors. *Journal of
##' Marketing Research, 18*(1), 39--50. \doi{10.2307/3151312}
##'
##' Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using
##' structural equation modeling: An alternative to coefficient alpha.
##' *Psychometrika, 74*(1), 155--167. \doi{10.1007/s11336-008-9099-3}
##'
##' McDonald, R. P. (1999). *Test theory: A unified treatment*. Mahwah, NJ:
##' Erlbaum.
##'
##' Raykov, T. (2001). Estimation of congeneric scale reliability using
##' covariance structure analysis with nonlinear constraints *British
##' Journal of Mathematical and Statistical Psychology, 54*(2), 315--323.
##' \doi{10.1348/000711001159582}
##'
##' Zumbo, B. D., Gadermann, A. M., & Zeisser, C. (2007). Ordinal versions of
##' coefficients alpha and theta for Likert rating scales.
##' *Journal of Modern Applied Statistical Methods, 6*(1), 21--29.
##' \doi{10.22237/jmasm/1177992180}
##'
##' Zumbo, B. D., & Kroc, E. (2019). A measurement is a choice and Stevens’
##' scales of measurement do not help make it: A response to Chalmers.
##' *Educational and Psychological Measurement, 79*(6), 1184--1197.
##' \doi{10.1177/0013164419844305}
##'
##'
##' @examples
##'
##' data(HolzingerSwineford1939)
##' HS9 <- HolzingerSwineford1939[ , c("x7","x8","x9")]
##' HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )
##' names(HSbinary) <- c("y7","y8","y9")
##' HS <- cbind(HolzingerSwineford1939, HSbinary)
##'
##' HS.model <- ' visual  =~ x1 + x2 + x3
##'               textual =~ x4 + x5 + x6
##'               speed   =~ y7 + y8 + y9 '
##'
##' fit <- cfa(HS.model, data = HS, ordered = c("y7","y8","y9"), std.lv = TRUE)
##'
##' ## works for factors with exclusively continuous OR categorical indicators
##' reliability(fit)
##'
##' ## reliability for ALL indicators only available when they are
##' ## all continuous or all categorical
##' reliability(fit, omit.factors = "speed", return.total = TRUE)
##'
##'
##' ## loop over visual indicators to calculate alpha if one indicator is removed
##' for (i in paste0("x", 1:3)) {
##'   cat("Drop x", i, ":\n")
##'   print(reliability(fit, omit.factors = c("textual","speed"),
##'                     omit.indicators = i, what = "alpha"))
##' }
##'
##'
##' ## works for multigroup models and for multilevel models (and both)
##' data(Demo.twolevel)
##' ## assign clusters to arbitrary groups
##' Demo.twolevel$g <- ifelse(Demo.twolevel$cluster %% 2L, "type1", "type2")
##' model2 <- ' group: type1
##'   level: within
##'     fac =~ y1 + L2*y2 + L3*y3
##'   level: between
##'     fac =~ y1 + L2*y2 + L3*y3
##'
##' group: type2
##'   level: within
##'     fac =~ y1 + L2*y2 + L3*y3
##'   level: between
##'     fac =~ y1 + L2*y2 + L3*y3
##' '
##' fit2 <- sem(model2, data = Demo.twolevel, cluster = "cluster", group = "g")
##' reliability(fit2, what = c("alpha","omega3"))
##'
##' @name reliability-deprecated
##' @usage
##' reliability(object, what = c("alpha", "omega", "omega2", "omega3", "ave"),
##'             return.total = FALSE, dropSingle = TRUE, omit.factors = character(0),
##'             omit.indicators = character(0), omit.imps = c("no.conv", "no.se"))
##' @seealso [semTools-deprecated()]
##' @keywords internal
NULL


##' @rdname semTools-deprecated
##' @section Reliability:
##' The original `reliability` function was suboptimally designed.
##' For example, AVE was returned, which is not a reliability index. Also,
##' alpha and several omega-type coefficients were returned, including the
##' original formula that was in appropriate for models with complex structure.
##' Some features could be controlled by the user for one but not both types of
##' index  For example, alpha for categorical indicators was returned on both
##' the observed and latent-response scales, but this was not an option for any
##' omega-type indices.  The omegas differed in terms of whether the observed or
##' model-implied covariance matrix was used in the denominator, but alpha was
##' only computed using the observed matrix.  These inconsistencies have been
##' resolved in the new [compRelSEM()] function, which returns only
##' one reliability index (per factor, optionally total score) according to the
##' user's requested features, for which there is much more flexibility.
##' Average variance extracted is now available in a dedicated [AVE()]
##' function.
##'
##' @export
reliability <- function(object,
                        what = c("alpha","omega","omega2","omega3","ave"),
                        return.total = FALSE, dropSingle = TRUE,
                        omit.factors = character(0),
                        omit.indicators = character(0),
                        omit.imps = c("no.conv","no.se")) {

  .Deprecated(msg = c("\nThe reliability() function was deprecated in 2022 and ",
                      "will cease to be included in future versions of semTools",
                      ". See help('semTools-deprecated) for details.",
                      "\n\nIt is replaced by the compRelSEM() function, which ",
                      "can estimate alpha and model-based reliability in an ",
                      "even wider variety of models and data types, with ",
                      "greater control in specifying the desired type of ",
                      "reliability coefficient (i.e., more explicitly choosing ",
                      "assumptions). \n\nThe average variance extracted should ",
                      "never have been included because it is not a reliability ",
                      "coefficient. It is now available from the AVE() function."))

  ngroups <- lavInspect(object, "ngroups") #TODO: adapt to multiple levels
  nLevels <- lavInspect(object, "nlevels")
  nblocks <- ngroups*nLevels #FIXME: always true?
  return.total <- rep(return.total, nblocks)
  group.label <- if (ngroups > 1L) lavInspect(object, "group.label") else NULL
  #FIXME? lavInspect(object, "level.labels")
  clus.label <- if (nLevels > 1L) c("within", lavInspect(object, "cluster")) else NULL
  if (nblocks > 1L) {
    block.label <- paste(rep(group.label, each = nLevels), clus.label,
                         sep = if (ngroups > 1L && nLevels > 1L) "_" else "")
  }

  ## check for categorical (determines what S will be)
  anyCategorical <- lavInspect(object, "categorical")
  if (anyCategorical && "alpha" %in% what) {
    what <- c(what, "alpha.ord")
    what <- unique(what) # in case it was already explicitly requested
  }
  ## categorical-model parameters
  threshold <- if (anyCategorical) getThreshold(object, omit.imps = omit.imps) else NULL
  latScales <- if (anyCategorical) getScales(object, omit.imps = omit.imps) else NULL
  ## all other relevant parameters in GLIST format (not flat, need block-level list)
  if (inherits(object, "lavaan")) {
    param <- lavInspect(object, "est")
    ve <- lavInspect(object, "cov.lv") # model-implied latent covariance matrix
    S <- object@h1$implied$cov # observed sample covariance matrix (already a list)
    if (anyCategorical && any(c("alpha","alpha.ord") %in% what)) {
      rawData <- try(lavInspect(object, "data"), silent = TRUE)
      if (inherits(rawData, "try-error"))
        stop('Error in lavInspect(fit, "data"); what="alpha" unavailable for ',
             'models fitted to summary statistics of categorial data.')
      if (nblocks == 1L) rawData <- list(rawData)
      S.as.con <- lapply(rawData, cov) # for actual "alpha", not "alpha.ord"
    }

    if (nblocks == 1L) {
      param <- list(param)
      ve <- list(ve)
    }

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    param <- object@coefList[[ useImps[1] ]] # first admissible as template
    coefList <- object@coefList[useImps]
    phiList <- object@phiList[useImps]
    if (anyCategorical) {
      dataList <- object@DataList[useImps]
      S.as.con <- vector("list", nblocks) # for group-list of pooled S
    }
    ## add block-level list per imputation?
    if (nblocks == 1L) {
      param <- list(param)
      for (i in 1:m) {
        coefList[[i]] <- list(coefList[[i]])
        phiList[[i]] <- list(phiList[[i]])
      }
      if (anyCategorical) { #FIXME: currently no categorical ML-SEMs
        #dataList[[i]] <- list(dataList[[i]])
        VV <- lavNames(object, type = "ov")
        impCovList <- lapply(dataList, function(DD) {
          dat <- do.call(cbind, sapply(DD[VV], as.numeric, simplify = FALSE))
          cov(dat)
        })
        S.as.con[[1]] <- Reduce("+", impCovList) / length(impCovList)
      }

    } else if (anyCategorical) { #FIXME: currently no categorical ML-SEMs
      ## multigroup models need separate data matrices per group
      G <- lavInspect(object, "group")

      for (g in seq_along(group.label)) {
        VV <- try(lavNames(object, type = "ov", group = group.label[g]),
                  silent = TRUE)
        if (inherits(VV, "try-error")) {
          VV <- lavNames(object, type = "ov", group = g)
        }
        impCovList <- lapply(dataList, function(DD) {
          RR <- DD[,G] == group.label[g]
          dat <- do.call(cbind, sapply(DD[RR, VV], as.numeric, simplify = FALSE))
          cov(dat)
        })
        S.as.con[[g]] <- Reduce("+", impCovList) / length(impCovList)
      }

    }
    S <- vector("list", nblocks) # pooled observed OR polychoric covariance matrix
    ve <- vector("list", nblocks)
    ## loop over blocks
    for (b in 1:nblocks) {

      ## param:  loop over GLIST elements
      for (mat in names(param[[b]])) {
        matList <- lapply(coefList, function(i) i[[b]][[mat]])
        param[[b]][[mat]] <- Reduce("+", matList) / length(matList)
      } # mat

      ## pooled observed OR polychoric covariance matrix
      covList <- lapply(object@h1List[useImps], function(i) i$implied$cov[[b]])
      S[[b]] <- Reduce("+", covList) / m

      ## pooled model-implied latent covariance matrix
      ve[[b]] <- Reduce("+", lapply(phiList, "[[", i = b) ) / m

    } # b

  }

  if (nblocks == 1L) {
    SigmaHat <- getMethod("fitted", class(object))(object)["cov"] # retain list format
  } else {
    SigmaHat <- sapply(getMethod("fitted", class(object))(object),
                       "[[", "cov", simplify = FALSE)
  }

  ly <- lapply(param, "[[", "lambda")
  te <- lapply(param, "[[", "theta")
  beta <- if ("beta" %in% names(param[[1]])) {
    lapply(param, "[[", "beta")
  } else NULL

	result <- list()
	warnTotal <- FALSE
	warnHigher <- character(0) # collect list of potential higher-order factors
	## loop over i blocks (groups/levels)
	for (i in 1:nblocks) {
	  ## extract factor and indicator names
	  allIndNames <- rownames(ly[[i]])
	  allFacNames <- colnames(ly[[i]])
	  myFacNames <- setdiff(allFacNames, omit.factors)
	  subLY <- ly[[i]][ , myFacNames, drop = FALSE] != 0
	  myIndNames <- rownames(subLY)[apply(subLY, MARGIN = 1L, FUN = any)]

	  ## distinguish between categorical, continuous, and latent indicators
	  nameArgs <- list(object = object)
	  if (nblocks > 1L) nameArgs$block <- i
	  ordNames <- do.call(lavNames, c(nameArgs, list(type = "ov.ord")))
	  numNames <- do.call(lavNames, c(nameArgs, list(type = "ov.num")))
	  if (anyCategorical) {
	    ## identify when the (sub)set of factors are all categorical
	    blockCat <- all(myIndNames %in% ordNames)
	    ## identify when the (sub)set of factors have mixed indicators, so no total
	    mix <- any(myIndNames %in% ordNames) && any(myIndNames %in% numNames)
	  } else {
	    blockCat <- FALSE
	    mix <- FALSE
	  }

	  if (mix && return.total[i]) {
	    return.total[i] <- FALSE
	    warnTotal <- TRUE
    }

	  ## identify POSSIBLE higher-order factors (that affect other latent vars)
	  latInds  <- do.call(lavNames, c(nameArgs, list(type = "lv.ind")))
	  higher <- if (length(latInds) == 0L) character(0) else {
	    allFacNames[apply(beta[[i]], MARGIN = 2, function(x) any(x != 0))]
	  }
	  ## keep track of factor indices to skip
	  idx.drop <- numeric(0)

	  ## relevant quantities
		common <- (apply(ly[[i]], 2, sum)^2) * diag(ve[[i]])
		truevar <- ly[[i]] %*% ve[[i]] %*% t(ly[[i]])
		## vectors to store results for each factor
		error <- rep(NA, length(common))
		alpha <- rep(NA, length(common))
		alpha.ord <- rep(NA, length(common))
		total <- rep(NA, length(common))
		omega1 <- omega2 <- omega3 <- rep(NA, length(common))
		impliedTotal <- rep(NA, length(common))
		avevar <- rep(NA, length(common))
		warnOmega <- FALSE
		## loop over j factors
		for (j in 1:length(common)) {
		  ## skip this factor?
		  if (allFacNames[j] %in% omit.factors) {
		    idx.drop <- c(idx.drop, j)
		    next
		  }

			index <- setdiff(which(ly[[i]][,j] != 0), # nonzero loadings
			                 which(allIndNames %in% omit.indicators))
			jIndNames <- allIndNames[index]

			## identify when this factor has mixed indicators, so no omegas
			jMix <- any(jIndNames %in% ordNames) && any(jIndNames %in% numNames)

			## check for ANY indicators (possibly skip purely higher-order factors)
			if (length(index) == 0L) {
			  idx.drop <- c(idx.drop, j)
			  next
			}
			## check for single indicators
			if (dropSingle && length(index) == 1L) {
			  idx.drop <- c(idx.drop, j)
			  next
			}
			## check for categorical (or mixed) indicators
			jCat <-      any(jIndNames %in% ordNames)
			warnOmega <- jCat && !all(jIndNames %in% ordNames)
			## check for latent indicators
			if (allFacNames[j] %in% higher && !(allFacNames[j] %in% omit.factors)) {
			  warnHigher <- c(warnHigher, allFacNames[j])
			}

			sigma <- S[[i]][index, index, drop = FALSE]
			faccontrib <- ly[[i]][,j, drop = FALSE] %*% ve[[i]][j,j, drop = FALSE] %*% t(ly[[i]][,j, drop = FALSE])
			truefac <- diag(faccontrib[index, index, drop = FALSE])
			trueitem <- diag(truevar[index, index, drop = FALSE])
			erritem <- diag(te[[i]][index, index, drop = FALSE])
			if (sum(abs(trueitem - truefac)) < 0.00001 & "ave" %in% what) {
				avevar[j] <- sum(trueitem) / sum(trueitem + erritem)
			}
			if (jCat) {
			  if ("alpha" %in% what) {
			    alpha[j] <- computeAlpha(S.as.con[[i]][index, index, drop = FALSE])
			  }
			  if ("alpha.ord" %in% what) {
			    alpha.ord[j] <- computeAlpha(sigma)
			  }
			  if ("omega" %in% what) {
			    omega1[j] <- omegaCat(truevar = faccontrib[index, index, drop = FALSE],
			                          threshold = threshold[[i]][jIndNames],
			                          scales = latScales[[i]][index],
			                          denom = faccontrib[index, index, drop = FALSE] + te[[i]][index, index, drop = FALSE])
			  }
			  if ("omega2" %in% what) {
			    omega2[j] <- omegaCat(truevar = faccontrib[index, index, drop = FALSE],
			                          threshold = threshold[[i]][jIndNames],
			                          scales = latScales[[i]][index],
			                          denom = SigmaHat[[i]][index, index, drop = FALSE])
			  }
			  if ("omega3" %in% what) {
			    omega3[j] <- omegaCat(truevar = faccontrib[index, index, drop = FALSE],
			                          threshold = threshold[[i]][jIndNames],
			                          scales = latScales[[i]][index],
			                          denom = sigma)
			  }

			} else {
			  alpha[j] <- computeAlpha(sigma)

			  commonfac <- sum(faccontrib[index, index, drop = FALSE])
			  error[j] <- sum(te[[i]][index, index, drop = FALSE])
			  impliedTotal[j] <- sum(SigmaHat[[i]][index, index, drop = FALSE])
			  total[j] <- sum(sigma)

			  omega1[j] <- commonfac / (commonfac + error[j])
				omega2[j] <- commonfac / impliedTotal[j]
				omega3[j] <- commonfac / total[j]
			}
			## end loop over j factors
		}

		if (return.total[i] & length(myFacNames) > 1L) {
		  if (blockCat) {
		    if ("alpha" %in% what) {
		      alpha <- c(alpha, computeAlpha(S.as.con[[i]]))
		    }
		    if ("alpha.ord" %in% what) {
		      alpha.ord <- c(alpha.ord, total = computeAlpha(S[[i]]))
		    }
		    if ("omega" %in% what) {
		      omega1 <- c(omega1, total = omegaCat(truevar = truevar,
		                                           threshold = threshold[[i]],
		                                           scales = latScales[[i]],
		                                           denom = truevar + te[[i]]))
		    }
		    if ("omega2" %in% what) {
		      omega2 <- c(omega2, total = omegaCat(truevar = truevar,
		                                           threshold = threshold[[i]],
		                                           scales = latScales[[i]],
		                                           denom = SigmaHat[[i]]))
		    }
		    if ("omega2" %in% what) {
  		    omega3 <- c(omega3, total = omegaCat(truevar = truevar,
  		                                         threshold = threshold[[i]],
  		                                         scales = latScales[[i]],
  		                                         denom = S[[i]]))
		    }

		  } else {
		    alpha <- c(alpha, total = computeAlpha(S[[i]]))
		    omega1 <- c(omega1, total = sum(truevar) / (sum(truevar) + sum(te[[i]])))
		    omega2 <- c(omega2, total = sum(truevar) / (sum(SigmaHat[[i]])))
		    omega3 <- c(omega3, total = sum(truevar) / (sum(S[[i]])))
		  }
		  avevar <- c(avevar,
		              total = sum(diag(truevar)) / sum((diag(truevar) + diag(te[[i]]))))
		}

		if (all(is.na(alpha.ord))) alpha.ord <- NULL
		result[[i]] <- rbind(alpha = if ("alpha" %in% what) alpha else NULL,
		                     alpha.ord = if ("alpha.ord" %in% what) alpha.ord else NULL,
		                     omega  = if ("omega"  %in% what) omega1 else NULL,
		                     omega2 = if ("omega2" %in% what) omega2 else NULL,
		                     omega3 = if ("omega3" %in% what) omega3 else NULL,
		                     avevar = if ("ave" %in% what) avevar else NULL)
		colnames(result[[i]])[1:length(allFacNames)] <- allFacNames
		if (return.total[i] & length(myFacNames) > 1L) {
		  colnames(result[[i]])[ ncol(result[[i]]) ] <- "total"
		}
		if (length(idx.drop)) {
		  result[[i]] <- result[[i]][ , -idx.drop, drop = FALSE]
		  ## reset indices for next block (could have different model/variables)
		  idx.drop <- numeric(0)
		}
		## end loop over blocks
	}

	warnCat <- sapply(result, function(x) any(c("alpha.ord","ave") %in% rownames(x)))
	if (any(warnCat)) {
	  alphaMessage <- paste0('Zumbo et al.`s (2007) "ordinal alpha" is calculated',
	                         ' in addition to the standard alpha, which treats ',
	                         'ordinal variables as numeric. See Chalmers (2018) ',
	                         'for a critique of "alpha.ord" and the response by ',
	                         'Zumbo & Kroc (2019).')
	  AVEmessage <- paste0('average variance extracted is calculated from ',
	                       'polychoric (polyserial) not Pearson correlations.')
	  both <- "alpha.ord" %in% what & "ave" %in% what
	  connectMessage <- if (both) ' Likewise, ' else ' the '
	  catMessage <- paste0("For constructs with categorical indicators, ",
	                       if ("alpha.ord" %in% what) alphaMessage else NULL,
	                       if (both) ' Likewise, ' else NULL,
	                       if ("ave" %in% what) AVEmessage else NULL)
	  if ("alpha.ord" %in% what || "ave" %in% what) message(catMessage, "\n")
	}
	if (length(warnHigher)) warning('Possible higher-order factors detected:\n',
	                                paste(unique(warnHigher), sep = ", "))
	if (warnTotal) {
	  message('Cannot return.total when model contains both continuous and ',
	          'binary/ordinal observed indicators. Use the ',
	          'omit.factors= argument to choose factors with only categorical ',
	          'indicators, if that is a composite of interest.\n')
	}
	if (warnOmega) {
	  message('Composite reliability (omega) cannot be computed for factors ',
	          'with mixed categorical and continuous indicators.')
	}

	## drop list structure?
	if (nblocks == 1L) {
		result <- result[[1]]
	} else names(result) <- block.label

	result
}



## ---------------
## reliabilityL2()
## (deprecated 10 May 2022)
## ---------------

##' Calculate the reliability values of a second-order factor
##'
##' Calculate the reliability values (coefficient omega) of a second-order
##' factor
##'
##' The first formula of the coefficient omega (in the
##' [reliability()]) will be mainly used in the calculation. The
##' model-implied covariance matrix of a second-order factor model can be
##' separated into three sources: the second-order common-factor variance,
##' the residual variance of the first-order common factors (i.e., not
##' accounted for by the second-order factor), and the measurement error of
##' observed indicators:
##'
##' \deqn{ \hat{\Sigma} = \Lambda \bold{B} \Phi_2 \bold{B}^{\prime}
##' \Lambda^{\prime} + \Lambda \Psi_{u} \Lambda^{\prime} + \Theta, }
##'
##' where \eqn{\hat{\Sigma}} is the model-implied covariance matrix,
##' \eqn{\Lambda} contains first-order factor loadings, \eqn{\bold{B}} contains
##' second-order factor loadings, \eqn{\Phi_2} is the covariance matrix of the
##' second-order factor(s), \eqn{\Psi_{u}} is the covariance matrix of residuals
##' from first-order factors, and \eqn{\Theta} is the covariance matrix of the
##' measurement errors from observed indicators. Thus, we can calculate the
##' proportion of variance of a composite score calculated from the observed
##' indicators (e.g., a total score or scale mean) that is attributable to the
##' second-order factor, i.e. coefficient omega at Level 1:
##'
##' \deqn{ \omega_{L1} = \frac{\bold{1}^{\prime} \Lambda \bold{B} \Phi_2
##' \bold{B}^{\prime} \Lambda^{\prime} \bold{1}}{\bold{1}^{\prime} \Lambda
##' \bold{B} \Phi_2 \bold{B} ^{\prime} \Lambda^{\prime} \bold{1} +
##' \bold{1}^{\prime} \Lambda \Psi_{u} \Lambda^{\prime} \bold{1} +
##' \bold{1}^{\prime} \Theta \bold{1}}, }
##'
##' where \eqn{\bold{1}} is the *k*-dimensional vector of 1 and *k* is
##' the number of observed variables.
##'
##' The model-implied covariance matrix among first-order factors (\eqn{\Phi_1})
##' can be calculated as:
##'
##' \deqn{ \Phi_1 = \bold{B} \Phi_2 \bold{B}^{\prime} + \Psi_{u}, }
##'
##' Thus, the proportion of variance among first-order common factors that is
##' attributable to the second-order factor (i.e., coefficient omega at Level 2)
##' can be calculated as:
##'
##' \deqn{ \omega_{L2} = \frac{\bold{1}_F^{\prime} \bold{B} \Phi_2
##' \bold{B}^{\prime} \bold{1}_F}{\bold{1}_F^{\prime} \bold{B} \Phi_2
##' \bold{B}^{\prime} \bold{1}_F + \bold{1}_F^{\prime} \Psi_{u} \bold{1}_F}, }
##'
##' where \eqn{\bold{1}_F} is the *F*-dimensional vector of 1 and *F*
##' is the number of first-order factors. This Level-2 omega can be interpreted
##' as an estimate of the reliability of a hypothetical composite calculated
##' from error-free observable variables representing the first-order common
##' factors. This might only be meaningful as a thought experiment.
##'
##' An additional thought experiment is possible: If the observed indicators
##' contained only the second-order common-factor variance and unsystematic
##' measurement error, then there would be no first-order common factors because
##' their unique variances would be excluded from the observed measures. An
##' estimate of this hypothetical composite reliability can be calculated as the
##' partial coefficient omega at Level 1, or the proportion of observed
##' variance explained by the second-order factor after partialling out the
##' uniqueness from the first-order factors:
##'
##' \deqn{ \omega_{L1} = \frac{\bold{1}^{\prime} \Lambda \bold{B} \Phi_2
##' \bold{B}^{\prime} \Lambda^{\prime} \bold{1}}{\bold{1}^{\prime} \Lambda
##' \bold{B} \Phi_2 \bold{B}^{\prime} \Lambda^{\prime} \bold{1} +
##' \bold{1}^{\prime} \Theta \bold{1}}, }
##'
##' Note that if the second-order factor has a direct factor loading on some
##' observed variables, the observed variables will be counted as first-order
##' factors, which might not be desirable.
##'
##'
##' @importFrom lavaan lavInspect
##'
##' @param object A [lavaan::lavaan-class] or [lavaan.mi::lavaan.mi-class] object,
##'   expected to contain a least one exogenous higher-order common factor.
##' @param secondFactor The name of a single second-order factor in the
##'   model fitted in `object`. The function must be called multiple
##'   times to estimate reliability for each higher-order factor.
##' @param omit.imps `character` vector specifying criteria for omitting
##'        imputations from pooled results.  Can include any of
##'        `c("no.conv", "no.se", "no.npd")`, the first 2 of which are the
##'        default setting, which excludes any imputations that did not
##'        converge or for which standard errors could not be computed.  The
##'        last option (`"no.npd"`) would exclude any imputations which
##'        yielded a nonpositive definite covariance matrix for observed or
##'        latent variables, which would include any "improper solutions" such
##'        as Heywood cases.  NPD solutions are not excluded by default because
##'        they are likely to occur due to sampling error, especially in small
##'        samples.  However, gross model misspecification could also cause
##'        NPD solutions, users can compare pooled results with and without
##'        this setting as a sensitivity analysis to see whether some
##'        imputations warrant further investigation.
##'
##' @return Reliability values at Levels 1 and 2 of the second-order factor, as
##'   well as the partial reliability value at Level 1
##'
##' @author Sunthud Pornprasertmanit (\email{psunthud@@gmail.com})
##'
##' @examples
##'
##' HS.model3 <- ' visual  =~ x1 + x2 + x3
##'                textual =~ x4 + x5 + x6
##'                speed   =~ x7 + x8 + x9
##'                higher =~ visual + textual + speed'
##'
##' fit6 <- cfa(HS.model3, data = HolzingerSwineford1939)
##' reliability(fit6) # Should provide a warning for the endogenous variables
##' reliabilityL2(fit6, "higher")
##'
##' @name reliabilityL2-deprecated
##' @usage
##' reliabilityL2(object, secondFactor, omit.imps = c("no.conv", "no.se"))
##' @seealso [semTools-deprecated()]
##' @keywords internal
NULL


##' @rdname semTools-deprecated
##' @section Higher-Order Reliability:
##' Originally, composite reliability of a single higher-order factor was
##' estimated in a separate function: `reliabilityL2`.  It is now available
##' for multiple higher-order factors in the same model, and from the same
##' [compRelSEM()] function that estimates reliability for first-order
##' factors, using the `higher=` argument to name higher-order factors in
##' the fitted model.
##'
##' @export
reliabilityL2 <- function(object, secondFactor,
                          omit.imps = c("no.conv","no.se")) {

  .Deprecated(msg = c("\nThe reliabilityL2() function was deprecated in 2022 and ",
                      "will cease to be included in future versions of semTools",
                      ". See help('semTools-deprecated) for details.",
                      "\n\nIt is replaced by the compRelSEM() function, which ",
                      "can estimate alpha and model-based reliability in an ",
                      "even wider variety of models and data types, with ",
                      "greater control in specifying the desired type of ",
                      "reliability coefficient (i.e., more explicitly choosing ",
                      "assumptions)."))

  secondFactor <- as.character(secondFactor)[1] # only one at a time

  ngroups <- lavInspect(object, "ngroups") #TODO: adapt to multiple levels
  nLevels <- lavInspect(object, "nlevels")
  nblocks <- ngroups*nLevels #FIXME: always true?
  group.label <- if (ngroups > 1L) lavInspect(object, "group.label") else NULL
  #FIXME? lavInspect(object, "level.labels")
  clus.label <- if (nLevels > 1L) c("within", lavInspect(object, "cluster")) else NULL
  if (nblocks > 1L) {
    block.label <- paste(rep(group.label, each = nLevels), clus.label,
                         sep = if (ngroups > 1L && nLevels > 1L) "_" else "")
  }

  ## parameters in GLIST format (not flat, need block-level list)
  if (inherits(object, "lavaan")) {
    param <- lavInspect(object, "est")
    ve <- lavInspect(object, "cov.lv") # model-implied latent covariance matrix
    S <- object@h1$implied$cov # observed sample covariance matrix (already a list)

    if (nblocks == 1L) {
      param <- list(param)
      ve <- list(ve)
    }

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    param <- object@coefList[[ useImps[1] ]] # first admissible as template
    coefList <- object@coefList[useImps]
    phiList <- object@phiList[useImps]
    ## add block-level list per imputation?
    if (nblocks == 1L) {
      param <- list(param)
      for (i in 1:m) {
        coefList[[i]] <- list(coefList[[i]])
        phiList[[i]] <- list(phiList[[i]])
      }
    }

    S <- vector("list", nblocks) # pooled observed covariance matrix
    ve <- vector("list", nblocks)
    ## loop over blocks
    for (b in 1:nblocks) {

      ## param:  loop over GLIST elements
      for (mat in names(param[[b]])) {
        matList <- lapply(coefList, function(i) i[[b]][[mat]])
        param[[b]][[mat]] <- Reduce("+", matList) / length(matList)
      } # mat

      ## pooled observed covariance matrix
      covList <- lapply(object@h1List[useImps], function(i) i$implied$cov[[b]])
      S[[b]] <- Reduce("+", covList) / m

      ## pooled model-implied latent covariance matrix
      ve[[b]] <- Reduce("+", lapply(phiList, "[[", i = b) ) / m

    } # b

  }

  if (nblocks == 1L) {
    SigmaHat <- getMethod("fitted", class(object))(object)["cov"] # retain list format
  } else {
    SigmaHat <- sapply(getMethod("fitted", class(object))(object),
                       "[[", "cov", simplify = FALSE)
  }

  ly <- lapply(param, "[[", "lambda")
  te <- lapply(param, "[[", "theta")
  ps <- lapply(param, "[[", "psi")
	be <- lapply(param, "[[", "beta")

	result <- list()
	for (i in 1:nblocks) {

		# Prepare for higher-order reliability
		l2var <- ve[[i]][secondFactor, secondFactor, drop = FALSE]
		l2load <- be[[1]][,secondFactor]
		indexl2 <- which(l2load != 0)
		commonl2 <- (sum(l2load)^2) * l2var
		errorl2 <- sum(ps[[i]][indexl2, indexl2, drop = FALSE])

		# Prepare for lower-order reliability
		indexl1 <- which(apply(ly[[i]][,indexl2, drop = FALSE], 1, function(x) sum(x != 0)) > 0)
		l1load <- ly[[i]][,indexl2] %*% as.matrix(be[[1]][indexl2, secondFactor, drop = FALSE])
		commonl1 <- (sum(l1load)^2) * l2var
		errorl1 <- sum(te[[i]][indexl1, indexl1, drop = FALSE])
		uniquel1 <- 0
		for (j in seq_along(indexl2)) {
			uniquel1 <- uniquel1 + (sum(ly[[i]][,indexl2[j]])^2) * ps[[i]][indexl2[j], indexl2[j], drop = FALSE]
		}

		# Adjustment for direct loading from L2 to observed variables
		if (any(ly[[i]][,secondFactor] != 0)) {
			indexind <- which(ly[[i]][,secondFactor] != 0)
			if (length(intersect(indexind, indexl1)) > 0)
			  stop("Direct and indirect loadings of higher-order factor to observed",
			       " variables are specified at the same time.")
			commonl2 <- sum(c(ly[[i]][,secondFactor], l2load))^2 * l2var
			errorl2 <- errorl2 + sum(te[[i]][indexind, indexind, drop = FALSE])
			commonl1 <- sum(c(ly[[i]][,secondFactor], l1load))^2 * l2var
			errorl1 <- errorl1 + sum(te[[i]][indexind, indexind, drop = FALSE])
		}

		# Calculate Reliability
		omegaL1 <- commonl1 / (commonl1 + uniquel1 + errorl1)
		omegaL2 <- commonl2 / (commonl2 + errorl2)
		partialOmegaL1 <- commonl1 / (commonl1 + errorl1)
		result[[i]] <- c(omegaL1 = omegaL1, omegaL2 = omegaL2, partialOmegaL1 = partialOmegaL1)
	}

	if (nblocks == 1L) {
	  result <- result[[1]]
	} else names(result) <- block.label

	result
}



## --------------
## maximalRelia()
## --------------

##' Calculate maximal reliability
##'
##' Calculate maximal reliability of a scale
##'
##' Given that a composite score (\eqn{W}) is a weighted sum of item scores:
##'
##' \deqn{ W = \bold{w}^\prime \bold{x} ,}
##'
##' where \eqn{\bold{x}} is a \eqn{k \times 1} vector of the scores of each
##' item, \eqn{\bold{w}} is a \eqn{k \times 1} weight vector of each item, and
##' \eqn{k} represents the number of items. Then, maximal reliability is
##' obtained by finding \eqn{\bold{w}} such that reliability attains its maximum
##' (Li, 1997; Raykov, 2012). Note that the reliability can be obtained by
##'
##' \deqn{ \rho = \frac{\bold{w}^\prime \bold{S}_T \bold{w}}{\bold{w}^\prime
##' \bold{S}_X \bold{w}}}
##'
##' where \eqn{\bold{S}_T} is the covariance matrix explained by true scores and
##' \eqn{\bold{S}_X} is the observed covariance matrix. Numerical method is used
##' to find \eqn{\bold{w}} in this function.
##'
##' For continuous items, \eqn{\bold{S}_T} can be calculated by
##'
##' \deqn{ \bold{S}_T = \Lambda \Psi \Lambda^\prime,}
##'
##' where \eqn{\Lambda} is the factor loading matrix and \eqn{\Psi} is the
##' covariance matrix among factors. \eqn{\bold{S}_X} is directly obtained by
##' covariance among items.
##'
##' For categorical items, Green and Yang's (2009) method is used for
##' calculating \eqn{\bold{S}_T} and \eqn{\bold{S}_X}. The element \eqn{i} and
##' \eqn{j} of \eqn{\bold{S}_T} can be calculated by
##'
##' \deqn{ \left[\bold{S}_T\right]_{ij} = \sum^{C_i - 1}_{c_i = 1} \sum^{C_j -
##' 1}_{c_j - 1} \Phi_2\left( \tau_{x_{c_i}}, \tau_{x_{c_j}}, \left[ \Lambda
##' \Psi \Lambda^\prime \right]_{ij} \right) - \sum^{C_i - 1}_{c_i = 1}
##' \Phi_1(\tau_{x_{c_i}}) \sum^{C_j - 1}_{c_j - 1} \Phi_1(\tau_{x_{c_j}}),}
##'
##' where \eqn{C_i} and \eqn{C_j} represents the number of thresholds in Items
##' \eqn{i} and \eqn{j}, \eqn{\tau_{x_{c_i}}} represents the threshold \eqn{c_i}
##' of Item \eqn{i}, \eqn{\tau_{x_{c_j}}} represents the threshold \eqn{c_i} of
##' Item \eqn{j}, \eqn{ \Phi_1(\tau_{x_{c_i}})} is the cumulative probability of
##' \eqn{\tau_{x_{c_i}}} given a univariate standard normal cumulative
##' distribution and \eqn{\Phi_2\left( \tau_{x_{c_i}}, \tau_{x_{c_j}}, \rho
##' \right)} is the joint cumulative probability of \eqn{\tau_{x_{c_i}}} and
##' \eqn{\tau_{x_{c_j}}} given a bivariate standard normal cumulative
##' distribution with a correlation of \eqn{\rho}
##'
##' Each element of \eqn{\bold{S}_X} can be calculated by
##'
##' \deqn{ \left[\bold{S}_T\right]_{ij} = \sum^{C_i - 1}_{c_i = 1} \sum^{C_j -
##' 1}_{c_j - 1} \Phi_2\left( \tau_{V_{c_i}}, \tau_{V_{c_j}}, \rho^*_{ij}
##' \right) - \sum^{C_i - 1}_{c_i = 1} \Phi_1(\tau_{V_{c_i}}) \sum^{C_j -
##' 1}_{c_j - 1} \Phi_1(\tau_{V_{c_j}}),}
##'
##' where \eqn{\rho^*_{ij}} is a polychoric correlation between Items \eqn{i}
##' and \eqn{j}.
##'
##'
##' @importFrom lavaan lavInspect lavNames
##'
##' @param object A [lavaan::lavaan-class] or [lavaan.mi::lavaan.mi-class] object,
##'   expected to contain only exogenous common factors (i.e., a CFA model).
##' @param omit.imps `character` vector specifying criteria for omitting
##'        imputations from pooled results.  Can include any of
##'        `c("no.conv", "no.se", "no.npd")`, the first 2 of which are the
##'        default setting, which excludes any imputations that did not
##'        converge or for which standard errors could not be computed.  The
##'        last option (`"no.npd"`) would exclude any imputations which
##'        yielded a nonpositive definite covariance matrix for observed or
##'        latent variables, which would include any "improper solutions" such
##'        as Heywood cases.  NPD solutions are not excluded by default because
##'        they are likely to occur due to sampling error, especially in small
##'        samples.  However, gross model misspecification could also cause
##'        NPD solutions, users can compare pooled results with and without
##'        this setting as a sensitivity analysis to see whether some
##'        imputations warrant further investigation.
##'
##' @return Maximal reliability values of each group. The maximal-reliability
##'   weights are also provided. Users may extracted the weighted by the
##'   `attr` function (see example below).
##'
##' @author Sunthud Pornprasertmanit (\email{psunthud@@gmail.com})
##'
##' @seealso [reliability()] for reliability of an unweighted
##'   composite score
##'
##' @references
##' Li, H. (1997). A unifying expression for the maximal reliability of a linear
##' composite. *Psychometrika, 62*(2), 245--249. \doi{10.1007/BF02295278}
##'
##' Raykov, T. (2012). Scale construction and development using structural
##' equation modeling. In R. H. Hoyle (Ed.), *Handbook of structural
##' equation modeling* (pp. 472--494). New York, NY: Guilford.
##'
##' @examples
##'
##' total <- 'f =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 '
##' fit <- cfa(total, data = HolzingerSwineford1939)
##' maximalRelia(fit)
##'
##' # Extract the weight
##' mr <- maximalRelia(fit)
##' attr(mr, "weight")
##'
##' @export
maximalRelia <- function(object, omit.imps = c("no.conv","no.se")) {
  ngroups <- lavInspect(object, "ngroups") #TODO: adapt to multiple levels
  nLevels <- lavInspect(object, "nlevels")
  nblocks <- ngroups*nLevels #FIXME: always true?
  group.label <- if (ngroups > 1L) lavInspect(object, "group.label") else NULL
  #FIXME? lavInspect(object, "level.labels")
  clus.label <- if (nLevels > 1L) c("within", lavInspect(object, "cluster")) else NULL
  if (nblocks > 1L) {
    block.label <- paste(rep(group.label, each = nLevels), clus.label,
                         sep = if (ngroups > 1L && nLevels > 1L) "_" else "")
  }

  ## parameters in GLIST format (not flat, need block-level list)
  if (inherits(object, "lavaan")) {
    param <- lavInspect(object, "est")
    ve <- lavInspect(object, "cov.lv") # model-implied latent covariance matrix
    S <- object@h1$implied$cov # observed sample covariance matrix (already a list)

    if (nblocks == 1L) {
      param <- list(param)
      ve <- list(ve)
    }

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    param <- object@coefList[[ useImps[1] ]] # first admissible as template
    coefList <- object@coefList[useImps]
    phiList <- object@phiList[useImps]
    ## add block-level list per imputation?
    if (nblocks == 1L) {
      param <- list(param)
      for (i in 1:m) {
        coefList[[i]] <- list(coefList[[i]])
        phiList[[i]] <- list(phiList[[i]])
      }
    }
    S <- vector("list", nblocks) # pooled observed covariance matrix
    ve <- vector("list", nblocks)
    ## loop over blocks
    for (b in 1:nblocks) {

      ## param:  loop over GLIST elements
      for (mat in names(param[[b]])) {
        matList <- lapply(coefList, function(i) i[[b]][[mat]])
        param[[b]][[mat]] <- Reduce("+", matList) / length(matList)
      } # mat

      ## pooled observed covariance matrix
      covList <- lapply(object@h1List[useImps], function(i) i$implied$cov[[b]])
      S[[b]] <- Reduce("+", covList) / m

      ## pooled model-implied latent covariance matrix
      ve[[b]] <- Reduce("+", lapply(phiList, "[[", i = b) ) / m

    } # b

  }

  if (nblocks == 1L) {
    SigmaHat <- getMethod("fitted", class(object))(object)["cov"] # retain list format
  } else {
    SigmaHat <- sapply(getMethod("fitted", class(object))(object),
                       "[[", "cov", simplify = FALSE)
  }

  ly <- lapply(param, "[[", "lambda")
  te <- lapply(param, "[[", "theta")

  categorical <- lavInspect(object, "categorical")
  threshold <- if (categorical) getThreshold(object, omit.imps = omit.imps) else NULL

  result <- list()
  for (i in 1:nblocks) {
    truevar <- ly[[i]] %*% ve[[i]] %*% t(ly[[i]])
    varnames <- colnames(truevar)
    if (categorical) {
      invstdvar <- 1 / sqrt(diag(SigmaHat[[i]]))
      polyr <- diag(invstdvar) %*% truevar %*% diag(invstdvar)
      nitem <- ncol(SigmaHat[[i]])
      result[[i]] <- calcMaximalReliaCat(polyr, threshold[[i]], S[[i]], nitem, varnames)
    } else {
      result[[i]] <- calcMaximalRelia(truevar, S[[i]], varnames)
    }
  }

  if (nblocks == 1L) {
    result <- result[[1]]
  } else names(result) <- block.label

  result
}



## ----------------
## Hidden Functions
## ----------------

computeAlpha <- function(S) {
  k <- nrow(S)
  k/(k - 1) * (1.0 - sum(diag(S)) / sum(S))
}

#' @importFrom stats cov2cor pnorm
omegaCat <- function(truevar, threshold, scales, denom) {
  ## must be in standardized latent scale
  R <- diag(scales) %*% truevar %*% diag(scales)

	## denom could be model-implied polychoric correlation assuming diagonal theta,
	##       model-implied polychoric correlation accounting for error covariances,
	##       or "observed" polychoric correlation matrix.
  ## If parameterization="theta", standardize the polychoric coVARIANCE matrix
	denom <- cov2cor(denom)

	nitem <- ncol(denom)
	## initialize sums of cumulative probabilities
	sumnum <- 0 # numerator
	addden <- 0 # denominator
	## loop over all pairs of items
	for (j in 1:nitem) {
  	for (jp in 1:nitem) {
  	  ## initialize sums of cumulative probabilities *per item*
  		sumprobn2 <- 0
  		addprobn2 <- 0
  		## for each pair of items, loop over all their thresholds
  		t1 <- threshold[[j]]  * scales[j] # on standardized latent scale
  		t2 <- threshold[[jp]] * scales[jp] #FIXME: subtract intercept (or marginal mean?)
  		for (c in 1:length(t1)) {
    		for (cp in 1:length(t2)) {
    			sumprobn2 <- sumprobn2 + p2(t1[c], t2[cp], R[j, jp])
    			addprobn2 <- addprobn2 + p2(t1[c], t2[cp], denom[j, jp])
    		}
  		}
  		sumprobn1 <- sum(pnorm(t1))
  		sumprobn1p <- sum(pnorm(t2))
  		sumnum <- sumnum + (sumprobn2 - sumprobn1 * sumprobn1p)
  		addden <- addden + (addprobn2 - sumprobn1 * sumprobn1p)
  	}
	}
	reliab <- sumnum / addden
	reliab
}


p2 <- function(t1, t2, r) {
	mnormt::pmnorm(c(t1, t2), c(0,0), matrix(c(1, r, r, 1), 2, 2))
}


# polycorLavaan <- function(object) {
# 	ngroups <- lavInspect(object, "ngroups")
# 	coef <- lavInspect(object, "est")
# 	targettaunames <- NULL
# 	if (ngroups == 1L) {
# 		targettaunames <- rownames(coef$tau)
# 	} else {
# 		targettaunames <- rownames(coef[[1]]$tau)
# 	}
# 	barpos <- sapply(strsplit(targettaunames, ""), function(x) which(x == "|"))
# 	varnames <- unique(apply(data.frame(targettaunames, barpos - 1), MARGIN = 1,
# 	                         FUN = function(x) substr(x[1], 1, x[2])))
# 	if (length(varnames))
# 	script <- ""
# 	for (i in 2:length(varnames)) {
# 		temp <- paste0(varnames[1:(i - 1)], collapse = " + ")
# 		temp <- paste0(varnames[i], "~~", temp, "\n")
# 		script <- paste(script, temp)
# 	}
# 	newobject <- refit(script, object)
# 	if (ngroups == 1L) {
# 		return(lavInspect(newobject, "est")$theta)
# 	}
# 	lapply(lavInspect(newobject, "est"), "[[", "theta")
# }

##' @importFrom lavaan lavInspect lavNames
getThreshold <- function(object, omit.imps = c("no.conv","no.se")) {
	ngroups <- lavInspect(object, "ngroups") #TODO: add nlevels when capable
	ordnames <- lavNames(object, "ov.ord")

	if (inherits(object, "lavaan")) {
	  EST <- lavInspect(object, "est")

	} else if (inherits(object, "lavaan.mi")) {
	  if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

	  useImps <- rep(TRUE, length(object@DataList))
	  if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
	  if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
	  if ("no.npd" %in% omit.imps) {
	    Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
	    Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
	    useImps <- useImps & !(Heywood.lv | Heywood.ov)
	  }
	  m <- sum(useImps)
	  if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
	  useImps <- which(useImps)

	  EST <- object@coefList[useImps]
	}

	if (ngroups == 1L) {
	  if (inherits(object, "lavaan")) {
	    thresholds <- EST$tau[,"threshold"]
	  } else if (inherits(object, "lavaan.mi")) {
	    tauList <- lapply(EST, function(x) x$tau[,"threshold"])
	    thresholds <- Reduce("+", tauList) / length(tauList)
	  }

	  result <- lapply(ordnames,
	                   function(nn) thresholds[grepl(paste0(nn, "\\|"), names(thresholds))])
	  names(result) <- ordnames
	  ## needs to be within a list when called above within block-loops
	  result <- list(result)

	} else {

	  thresholds <- vector("list", ngroups)
	  for (g in 1:ngroups) {
	    if (inherits(object, "lavaan")) {
	      thresholds[[g]] <- EST[[g]]$tau[,"threshold"]
	    } else if (inherits(object, "lavaan.mi")) {
	      tauList <- lapply(EST, function(x) x[[g]]$tau[,"threshold"])
	      thresholds[[g]] <- Reduce("+", tauList) / length(tauList)
	    }

	  }

	  result <- list()
		group.label <- lavInspect(object, "group.label")

		for (g in 1:ngroups) {
		  result[[ group.label[g] ]] <- lapply(ordnames, function(nn) {
		    thresholds[[g]][ grepl(paste0(nn, "\\|"), names(thresholds[[g]])) ]
		  })
		  names(result[[ group.label[g] ]]) <- ordnames
		}

	}

	return(result)
}

##' @importFrom lavaan lavInspect lavNames
getScales <- function(object, omit.imps = c("no.conv","no.se")) {
  ngroups <- lavInspect(object, "ngroups") #TODO: add nlevels when capable
  ordnames <- lavNames(object, "ov.ord") #TODO: use to allow mix of cat/con vars

  if (inherits(object, "lavaan")) {
    EST <- lavInspect(object, "est")

  } else if (inherits(object, "lavaan.mi")) {
    if (!"package:lavaan.mi" %in% search()) attachNamespace("lavaan.mi")

    useImps <- rep(TRUE, length(object@DataList))
    if ("no.conv" %in% omit.imps) useImps <- sapply(object@convergence, "[[", i = "converged")
    if ("no.se" %in% omit.imps) useImps <- useImps & sapply(object@convergence, "[[", i = "SE")
    if ("no.npd" %in% omit.imps) {
      Heywood.lv <- sapply(object@convergence, "[[", i = "Heywood.lv")
      Heywood.ov <- sapply(object@convergence, "[[", i = "Heywood.ov")
      useImps <- useImps & !(Heywood.lv | Heywood.ov)
    }
    m <- sum(useImps)
    if (m == 0L) stop('No imputations meet "omit.imps" criteria.')
    useImps <- which(useImps)

    EST <- object@coefList[useImps]
  }

  if (ngroups == 1L) {
    if (inherits(object, "lavaan")) {
      result <- list(EST$delta[,"scales"])
    } else if (inherits(object, "lavaan.mi")) {
      scales <- lapply(EST, function(x) x$delta[,"scales"])
      result <- list(Reduce("+", scales) / length(scales))
    }

  } else {

    result <- vector("list", ngroups)

    for (g in 1:ngroups) {
      if (inherits(object, "lavaan")) {
        result[[g]] <- EST[[g]]$delta[,"scales"]
      } else if (inherits(object, "lavaan.mi")) {
        scales <- lapply(EST, function(x) x[[g]]$delta[,"scales"])
        result[[g]] <- Reduce("+", scales) / length(scales)
      }
    }

  }

  return(result)
}

invGeneralRelia <- function(w, truevar, totalvar) {
	1 - (t(w) %*% truevar %*% w) / (t(w) %*% totalvar %*% w)
}

#' @importFrom stats pnorm
invGeneralReliaCat <- function(w, polyr, threshold, denom, nitem) {
	# denom could be polychoric correlation, model-implied correlation, or model-implied without error correlation
	upper <- matrix(NA, nitem, nitem)
	lower <- matrix(NA, nitem, nitem)
	for (j in 1:nitem) {
  	for (jp in 1:nitem) {
  		sumprobn2 <- 0
  		addprobn2 <- 0
  		t1 <- threshold[[j]]
  		t2 <- threshold[[jp]]
  		for (c in 1:length(t1)) {
  		for (cp in 1:length(t2)) {
  			sumprobn2 <- sumprobn2 + p2(t1[c], t2[cp], polyr[j, jp])
  			addprobn2 <- addprobn2 + p2(t1[c], t2[cp], denom[j, jp])
  		}
  		}
  		sumprobn1 <- sum(pnorm(t1))
  		sumprobn1p <- sum(pnorm(t2))
  		upper[j, jp] <- (sumprobn2 - sumprobn1 * sumprobn1p)
  		lower[j, jp] <- (addprobn2 - sumprobn1 * sumprobn1p)
  	}
	}
	1 - (t(w) %*% upper %*% w) / (t(w) %*% lower %*% w)
}

#' @importFrom stats nlminb
calcMaximalRelia <- function(truevar, totalvar, varnames) {
	start <- rep(1, nrow(truevar))
	out <- nlminb(start, invGeneralRelia, truevar = truevar, totalvar = totalvar)
	if (out$convergence != 0) stop("The numerical method for finding the maximal",
	                               " reliability did not converge.")
	result <- 1 - out$objective
	weight <- out$par / mean(out$par)
	names(weight) <- varnames
	attr(result, "weight") <- weight
	result
}

#' @importFrom stats nlminb
calcMaximalReliaCat <- function(polyr, threshold, denom, nitem, varnames) {
	start <- rep(1, nrow(polyr))
	out <- nlminb(start, invGeneralReliaCat, polyr = polyr, threshold = threshold,
	              denom = denom, nitem = nitem)
	if (out$convergence != 0) stop("The numerical method for finding the maximal",
	                               " reliability did not converge.")
	result <- 1 - out$objective
	weight <- out$par / mean(out$par)
	names(weight) <- varnames
	attr(result, "weight") <- weight
	result
}