Nothing
R/multilevel.omega.R
In misty: Miscellaneous Functions 'T. Yanagida'
Defines functions
multilevel.omega
Documented in
multilevel.omega
#' Multilevel Composite Reliability
#'
#' This function computes point estimate and Monte Carlo confidence interval for
#' the multilevel composite reliability defined by Lai (2021) for a (1) within-cluster
#' construct, (2) shared cluster-level construct, and (3) individual and configural
#' cluster construct by calling the \code{cfa} function in the R package \pkg{lavaan}.
#' By default, the function prints level-specific multilevel composite reliability
#' indices for an individual and configural cluster construct with 95% Monte Carlo
#' confidence interval based on Huber-White standard errors.
#'
#' @param data a data frame. Multilevel confirmatory factor analysis
#' based on a measurement model with one factor at the Within
#' level and one factor at the Between level comprising all
#' variables in the data frame is conducted. Note that the
#' cluster variable specified in \code{cluster} is excluded
#' from \code{data} when specifying the argument \code{cluster}
#' using the variable name of the cluster variable.
#' @param ... an expression indicating the variable names in \code{data},
#' e.g., \code{multilevel.omega(dat, x1, x2, x3, cluster = "cluster)}.
#' Note that the operators \code{+}, \code{-},
#' \code{~}, \code{:}, \code{::}, and \code{!} can also be
#' used to select variables, see 'Details' in the
#' \code{\link{df.subset}} function.
#' @param cluster either a character string indicating the variable name of
#' the cluster variable in \code{data}, or a vector representing
#' the nested grouping structure (i.e., group or cluster variable).
#' @param rescov a character vector or a list of character vectors for specifying
#' residual covariances at the Within level, e.g. \code{rescov = c("x1", "x2")}
#' for specifying a residual covariance between indicators \code{x1}
#' and \code{x2} at the Within level or \code{rescov = list(c("x1", "x2"), c("x3", "x4"))}
#' for specifying residual covariances between indicators \code{x1}
#' and \code{x2}, and indicators \code{x3} and \code{x4} at
#' the Within level. Note that residual covariances at the
#' Between level cannot be specified using this function.
#' @param const a character string indicating the type of construct(s), i.e.,
#' \code{"within"} for within-cluster constructs, \code{"shared"}
#' for shared cluster-level constructs, and \code{"config"}
#' (default) for individual and configural constructs.
#' @param fix.resid a character vector for specifying residual variances to be
#' fixed at 0 at the Between level, e.g., \code{fix.resid = c("x1", "x3")}
#' to fix residual variances of indicators \code{x1} and \code{x2}
#' at the Between level at 0. Note that it is also possible
#' to specify \code{fix.resid = "all"} which fixes all residual
#' variances at the Between level at 0 in line with the strong
#' factorial measurement invariance assumption across cluster.
#' @param se a character string indicating the standard errors used for
#' computing Monte Carlo confidence intervals, i.e.,
#' \code{"none"} for no standard errors, \code{"standard"} for
#' conventional standard error based on inverting the expected
#' observed or first.order information matrix, and
#' \code{"robust.huber.white"} for the 'MLR' (aka pseudo ML,
#' Huber-White) approach. Note that \code{se = "none"} saves
#' computation time by not computing standard errors and
#' Monte Carlo confidence intervals.
#' @param optim.method a character string indicating the optimizer, i.e., \code{"nlminb"}
#' (default) for the unconstrained and bounds-constrained
#' quasi-Newton method optimizer and \code{"em"} for the
#' Expectation Maximization (EM) algorithm.
#' @param missing a character string indicating how to deal with missing data,
#' i.e., \code{"listwise"} for listwise deletion or \code{"fiml"}
#' (default) for full information maximum likelihood (FIML)
#' method.
#' @param nrep an integer value indicating the number of Monte Carlo
#' repetitions for computing confidence intervals.
#' @param seed a numeric value specifying the seed of the random number
#' generator for computing the Monte Carlo confidence interval.
#' @param conf.level a numeric value between 0 and 1 indicating the confidence
#' level of the interval.
#' @param print a character vector indicating which results to show, i.e.
#' \code{"all"} (default) for all results \code{"omega"} for
#' the composite reliability omega and \code{"item"} for item
#' statistics.
#' @param digits an integer value indicating the number of decimal places
#' to be used for displaying mean, standard deviation, minimum,
#' maximum, skewness, and kurtosis.
#' @param r.digits an integer value indicating the number of decimal places
#' to be used for displaying multilevel coefficient omega,
#' ICC(1), and standardized factor loadings.
#' @param as.na a numeric vector indicating user-defined missing values,
#' i.e. these values are converted to \code{NA} before conducting
#' the analysis. Note that \code{as.na()} function is only
#' applied to \code{data} but not to \code{cluster}.
#' @param write a character string naming a file for writing the output into
#' either a text file with file extension \code{".txt"} (e.g.,
#' \code{"Output.txt"}) or Excel file with file extension
#' \code{".xlsx"} (e.g., \code{"Output.xlsx"}). If the file
#' name does not contain any file extension, an Excel file will
#' be written.
#' @param append logical: if \code{TRUE} (default), output will be appended
#' to an existing text file with extension \code{.txt} specified
#' in \code{write}, if \code{FALSE} existing text file will be
#' overwritten.
#' @param check logical: if \code{TRUE} (default), argument specification,
#' convergence and model identification is checked.
#' @param output logical: if \code{TRUE} (default), output is shown.
#'
#' @details
#' Geldhof et al. (2014) introduced a multilevel confirmatory factor analysis
#' approach for estimating score reliability at each level. This approach uses
#' the estimated factor loadings, factor variances, and residual variances from both
#' the within-level and between-level model to compute the within-level reliability
#' \eqn{\tilde{\omega}^w} and the between-level reliability \eqn{\tilde{\omega}^b}
#' based on McDonald's \eqn{\omega}. However, Lai (2021) pointed out two limitations
#' of this approach. First, it does not account for different conceptualizations
#' of constructs in multilevel data, where constructs may have different meanings
#' (Stapleton et al., 2016). Second, the between-level reliability \eqn{\tilde{\omega}^b}
#' in Geldhof et al. (2014) is a measure of reliability of the latent cluster means,
#' which ignores the sampling error in the observed cluster means (see Lüdtke et
#' al., 2011). Consequently, \eqn{\tilde{\omega}^b} substantially overestimates
#' the true reliability of observed between-level composite scores. To address
#' these limitations, Lai (2021) proposed multilevel reliability coefficients
#' tailored to different types of constructs to estimate the reliability of the
#' observed scores at each level: (1) within-cluster construct, (2) shared
#' cluster-level construct, and (3) individual and configural constructs.
#'
#' \describe{
#' \item{\strong{Within-Cluster Construct}}{The latent variable representing the
#' within-cluster construct is only meaningful at the within level. For example,
#' students' popularity among their peers within the same classroom based on
#' sociometric ratings. In this case any variation across classrooms should be
#' irrelevant to the construct. Thus, a saturated between-level model is specified.
#' The reliability for a within-level composite, \eqn{Z^w_{ij}}, measuring
#' a within-cluster construct is
#'
#' \deqn{\omega^w = \frac{ (\sum^p_{k = 1}\lambda^w_k)^2\phi^w} {(\sum^p_{k = 1}\lambda^w_k)^2\phi^w + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' Note that the formula is the same as the formula for the individual construct
#' of an individual and configural construct, except that the factor loadings are
#' specific to the within level as the between-level model is saturated.}
#'
#' \item{\strong{Shared Cluster-Level Construct}}{The latent variable representing
#' the shared construct is only meaningful at the between level. For example,
#' safety of a neighborhood, effectiveness of a teacher, or organizational
#' climate are measured by multiple informants at the within level. In these cases
#' any variation within clusters should be irrelevant to the construct. Thus, a
#' saturated within-level model is specified. The reliability for a
#' between-level composite, \eqn{Z^b_{j}}, measuring a shared construct is
#'
#' \deqn{\omega^b = \frac{(\sum^p_{k = 1}\lambda^b_k)^2\phi^b} {(\sum^p_{k = 1}\lambda^b_k)^2\phi^b + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1}} + \mathbf{1}^{\prime}\mathbf{\sum}^w\mathbf{1}/\tilde{n}}
#'
#' where \eqn{\mathbf{\sum}^w} is the variance covariance matrix of the saturated
#' within-level model and \eqn{\tilde{n}} is the harmonic mean of cluster sizes
#' \eqn{\tilde{n} = \frac{1}{\sum^p_{k = 1}\frac{1}{n_j}}}.}
#'
#' \item{\strong{Individual and Configural Construct}}{The individual construct is
#' defined at the within level (e.g., individual student achievement), whereas the
#' configural construct consists of cluster averages of individual constructs (e.g.,
#' mean student achievement of a classroom). Note that the factor loadings need to
#' be constrained to be equal across levels, i.e., \eqn{\mathbf{\lambda}^w = \mathbf{\lambda}^b}
#' so that the latent variables at the within level, \eqn{\eta^w}, and between level,
#' \eqn{\eta^b}, are on the same metric.
#'
#' \describe{
#' \item{\emph{Individual Construct}}{The reliability for a within-level
#' composite, \eqn{Z^w_{ij}}, measuring an individual construct is
#'
#' \deqn{\omega^w = \frac{ (\sum^p_{k = 1}\lambda^k)^2\phi^w} {(\sum^p_{k = 1}\lambda^k)^2\phi^w + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' }
#'
#' \item{\emph{Configural Construct}}{The reliability for a between-level
#' composite, \eqn{Z^b_{j}}, measuring a configural construct is
#'
#' \deqn{\omega^b = \frac{(\sum^p_{k = 1}\lambda_k)^2\phi^b}{(\sum^p_{k = 1}\lambda_k)^2(\phi^b + \phi^w / \tilde{n}) + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1} + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}/\tilde{n}}}}
#'
#' \item{\emph{Composite Measuring an Individual and Configural Construct}}{
#' The reliability for an overall composite, \eqn{Z_{ij}}, measuring an individual
#' construct capturing the population variance of both the within-level and
#' the between-level components of the true score and the errors is
#'
#' \deqn{\omega^{2l} = \frac{(\sum^p_{k = 1}\lambda_k)^2(\phi^w + \phi^b)}{(\sum^p_{k = 1}\lambda_k)^2(\phi^w + \phi^b) + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1} + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' Note that \eqn{\omega^{2l}} is simply the ratio of true variance to total
#' variance of the observed scores. In practice, \eqn{\omega^{2l}} is not
#' the most interesting coefficient because researchers are mainly interested
#' in distinguishing within-and between effects (Castro-Alvarez et al., 2026).}
#'
#' }
#' }
#' }
#' Note that the multilevel factor model assumes \eqn{p} items measuring one latent
#' construct at the within level (within-cluster construct), between level (shared
#' cluster-level construct), or within- and between level (individual and configural
#' construct). Local independence is assumed, \eqn{\Theta^w_j = diag[\theta^w_{11} ... \theta^w_{pp}]}
#' and \eqn{\Theta^b_j = diag[\theta^b_{11} ... \theta^b_{pp}]}, but composite reliability
#' can be computed when this assumption is violated by specifying the \code{rescov}
#' argument. The model assumes equal factor loadings across clusters,
#' \eqn{\mathbf{\lambda}^w_j = \mathbf{\lambda}^w}, and homogeneity of error
#' covariances across cluster, \eqn{\mathbf{\Theta}^w_j = \mathbf{\Theta}^w}.
#' In case of an individual and configural construct, the model also assumes cross-level
#' measurement invariance, i.e., \eqn{\mathbf{\lambda}^w = \mathbf{\lambda}^b}}.
#'
#' @author
#' Takuya Yanagida \email{takuya.yanagida@@univie.ac.at}
#'
#' @seealso
#' \code{\link{multilevel.alpha}}, \code{\link{item.omega}}, \code{\link{multilevel.cfa}},
#' \code{\link{multilevel.fit}}, \code{\link{multilevel.invar}}, \code{\link{multilevel.cor}},
#' \code{\link{multilevel.descript}}, \code{\link{write.result}}
#'
#' @references
#' Castro-Alvarez, S., Bringmann, L. F., Back, J., & Liu, S. (2026). The many
#' reliabilities of psychological dynamics: An overview of statistical approaches
#' to estimate the internal consistency reliability of intensive longitudinal data.
#' \emph{Psychological Methods, 31}(2), 281–296. https://doi.org/10.1037/met0000778
#'
#' Geldhof, G. J., Preacher, K. J., & Zyphur, M. J. (2014). Reliability
#' estimation in a multilevel confirmatory factor analysis framework. \emph{Psychological
#' Methods, 19}, 72-91. http://dx.doi.org/10.1037/a0032138
#'
#' Lai, M. H. C. (2021). Composite reliability of multilevel data: It’s about
#' observed scores and construct meanings. \emph{Psychological Methods, 26}(1),
#' 90–102. https://doi.org/10.1037/met0000287
#'
#' Lüdtke, O., Marsh, H. W., Robitzsch, A., & Trautwein, U. (2011). A 2 2
#' taxonomy of multilevel latent contextual models: Accuracy bias tradeoffs
#' in full and partial error correction models. \emph{Psychological Methods, 16},
#' 444-467. http://dx.doi.org/10.1037/a0024376
#'
#' Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling.
#' \emph{Journal of Statistical Software, 48}, 1-36. https://doi.org/10.18637/jss.v048.i02
#'
#' Stapleton, L. M., Yang, J. S., & Hancock, G. R. (2016). Construct meaning in
#' multilevel settings. \emph{Journal of Educational and Behavioral Statistics, 41}(5),
#' 481–520. https://doi.org/10.3102/1076998616646200
#'
#' Venables, W. N., Ripley, B. D. (2002).\emph{Modern Applied Statistics with S} (4th ed.).
#' Springer. https://www.stats.ox.ac.uk/pub/MASS4/.
#'
#' @return
#' Returns an object of class \code{misty.object}, which is a list with following
#' entries:
#' \item{\code{call}}{function call}
#' \item{\code{type}}{type of analysis}
#' \item{\code{data}}{data frame specified in \code{data} including the group variable
#' specified in \code{cluster}}
#' \item{\code{args}}{specification of function arguments}
#' \item{\code{model}}{specified model}
#' \item{\code{model.fit}}{fitted lavaan object (\code{mod.fit})}
#' \item{\code{check}}{results of the convergence and model identification check}
#' \item{\code{result}}{list with result tables, i.e., \code{omega} for the coefficient
#' omega including Monte Carlo confidence interval and
#' \code{itemstat} for descriptive statistics}
#'
#' @note
#' The function uses the functions \code{lavInspect}, \code{lavTech}, and \code{lavNames},
#' provided in the R package \pkg{lavaan} by Yves Rosseel (2012). The internal function
#' \code{.internal.mvrnorm} is a copy of the \code{mvrnorm} function in the package
#' \pkg{MASS} by Venables and Ripley (2002).
#'
#' @export
#'
#' @examples
#' \dontrun{
#' # Load data set "Demo.twolevel" in the lavaan package
#' data("Demo.twolevel", package = "lavaan")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Cluster Variable Specification
#'
#' # Example 1a: Specification using the argument '...'
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster")
#'
#' # Example 1b: Alternative specification with cluster variable 'cluster' in 'data'
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4", "cluster")], cluster = "cluster")
#'
#' # Example 1b: Alternative specification with cluster variable 'cluster' not in 'data'
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")], cluster = Demo.twolevel$cluster)
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Type of Construct
#'
#' # Example 2a: Within-Cluster Construct
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")],
#' cluster = Demo.twolevel$cluster, const = "within")
#'
#' # Example 2b: Shared Cluster-Level Construct
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "shared")
#'
#' # Example 2c: Configural Construct
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Residual (Co)Variances
#'
#' # Example 3a: Residual covariance between "y4" and "y5" at the Within level
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config",
#' rescov = c("y3", "y4"))
#'
#' # Example 3b: Residual variances of 'y1' at the Between level fixed at 0
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config",
#' fix.resid = c("y1", "y2"))
#'
#' #————————————————————————————————————————————————————————————————————————————
# # Arguments 'se' and 'print'
#'
#' # Example 4a: No confidence intervals to speed up computation
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster", se = "none")
#'
#' # Example 4b: Omega and item statistics
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster", print = "all")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Write Results
#'
#' # Example 4a: Write results into a text file
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")],
#' cluster = Demo.twolevel$cluster, write = "Multilevel_Omega.txt")
#'
#' # Example 4b: Write results into a Excel file
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster",
#' write = "Multilevel_Omega.xlsx")
#' }
multilevel.omega <- function(data, ..., cluster, rescov = NULL,
const = c("within", "shared", "config"),
fix.resid = NULL, se = c("none", "standard", "robust.huber.white"),
optim.method = c("nlminb", "em"), missing = c("listwise", "fiml"),
nrep = 100000, seed = NULL, conf.level = 0.95,
print = c("all", "omega", "item"), digits = 2, r.digits = 3,
as.na = NULL, write = NULL, append = TRUE, check = TRUE,
output = TRUE) {
#_____________________________________________________________________________
#
# Initial Check --------------------------------------------------------------
# Check if input 'data' is missing or NULL
if (isTRUE(missing(data) || is.null(data))) { stop("Please specify a data frame for the argument 'data'", call. = FALSE) }
# Check input 'cluster'
if (isTRUE(missing(cluster) ||is.null(cluster))) { stop("Please specify a variable name or vector representing the grouping structure for the argument 'cluster'.", call. = FALSE) }
#_____________________________________________________________________________
#
# Data -----------------------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Using the Argument '...' ####
if (isTRUE(!missing(...))) {
# Extract data
x <- as.data.frame(data[, .var.names(data = data, ..., cluster = cluster), drop = FALSE])
# Extract cluster variable and convert tibble into data frame or vector
cluster <- data[, cluster] |> (\(y) if (isTRUE("tbl" %in% substr(class(y), 1L, 3L))) { unname(unlist(y)) } else { return(y) })()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Without Using the Argument '...' ####
} else {
# Data frame
x <- as.data.frame(data)
# Data and cluster
var.group <- .var.group(data = x, cluster = cluster)
# Data
if (isTRUE(!is.null(var.group$data))) { x <- var.group$data }
# Cluster variable
if (isTRUE(!is.null(var.group$cluster))) { cluster <- var.group$cluster }
}
#_____________________________________________________________________________
#
# Input Check ----------------------------------------------------------------
# Check inputs
.check.input(logical = c("append", "output"),
s.character = list(const = c("within", "shared", "config"), se = c("none", "standard", "robust.huber.white"), optim.method = c("nlminb", "em"), missing = c("listwise", "fiml")),
m.character = list(print = c("all", "omega", "item")),
args = c("conf.level", "digits", "r.digits", "nrep", "seed", "write2"), envir = environment(), input.check = check)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Additional Checks
if (isTRUE(check)) {
#—————————————————————————————————————— #
### Check Input 'rescov' ####
if (isTRUE(!is.null(rescov))) {
# Two variables for each residual covariance
if (isTRUE(is.list(rescov) && any(sapply(rescov, length) != 2L))) { stop("Please specify a list of character vectors for the argument 'rescov', where each element has two variable names", call. = FALSE)
} else { if (isTRUE(length(rescov) != 2L)) { stop("Please specify a character vector with two variable names for the argument 'rescov'", call. = FALSE) } }
# Variable in 'data'
(!unique(unlist(rescov)) %in% colnames(x)) |> (\(y) if (isTRUE(any(y))) { stop(paste0("Variables specified in the argument 'rescov' were not found in 'data': ", paste(unique(unlist(rescov))[y], collapse = ", ")), call. = FALSE) })()
}
#—————————————————————————————————————— #
### Check Input 'fix.resid' ####
(!unique(fix.resid) %in% colnames(x)) |> (\(y) if (isTRUE(any(y) && all(fix.resid != "all"))) { stop(paste0("Variables specified in the argument 'fix.resid' were not found in 'data': ", paste(fix.resid[y], collapse = ", ")), call. = FALSE) })()
}
#_____________________________________________________________________________
#
# Data and Arguments ---------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Manifest Variables ####
var <- colnames(x)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Data Frame with Cluster Variable ####
x <- data.frame(x, .cluster = cluster, row.names = NULL)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing Data on the Cluster Variable ####
if (isTRUE(any(is.na(x$.cluster)))) {
warning(paste0("Data contains missing values on the cluster variable, number of cases removed from the analysis: ", sum(is.na(x$.cluster))), call. = FALSE)
x <- x[!is.na(x$.cluster), ]
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Convert User-Missing Values into NA ####
if (isTRUE(!is.null(as.na))) { x[, var] <- .as.na(x[, var], na = as.na) }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Type of Construct(s) ####
if (isTRUE(all(c("within", "shared", "config") %in% const))) { const <- "config" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Residual Variances Fixed at 0 ####
if (isTRUE(fix.resid == "all")) { fix.resid <- var }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Standard Error ####
if (isTRUE(all(c("none", "standard", "robust.huber.white") %in% se))) { se <- "robust.huber.white" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Optimizer ####
if (isTRUE(all(c("nlminb", "em") %in% optim.method))) { optim.method <- "nlminb" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing ####
# Complete data
if (isTRUE(all(!is.na(x[, var])))) {
missing <- "listwise"
# Data with missing values
} else {
if (isTRUE(all(c("listwise", "fiml") %in% missing))) {
missing <- "fiml"
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing Data on All Variable ####
(misty::na.prop(x[, var], append = FALSE) == 1L) |> (\(y) if (isTRUE(any(y) && missing == "fiml")) { warning(paste0("Data contains cases with missing values on all variables, number of cases removed from the analysis: ", sum(y)), call. = FALSE) })()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Print ####
if (isTRUE(all(c("all", "omega", "item") %in% print))) {
print <- "omega"
} else if (isTRUE(all(print == "all"))) {
print <- c("omega", "item")
}
#_____________________________________________________________________________
#
# Main Function --------------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Model Estimation ####
model.fit <- tryCatch(suppressWarnings(misty::multilevel.cfa(x, cluster = ".cluster", model = NULL, rescov = rescov,
model.w = NULL, model.b = NULL, rescov.w = NULL, rescov.b = NULL,
const = const, fix.resid = fix.resid, ident = "var", ls.fit = FALSE,
estimator = "ML", optim.method = optim.method, test = "none", se = se,
missing = missing, output = FALSE, check = FALSE)),
error = function(y) {
if (isTRUE(missing == "fiml")) {
stop("There was an estimation problem in lavaan, switching to missing = \"listwise\" might solve the problem.", call. = FALSE)
} else {
stop("There was an estimation problem in lavaan, model could not be estimated.", call. = FALSE)
}})
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Convergence and Model Identification Checks ####
if (isTRUE(check)) {
#—————————————————————————————————————— #
### Model Convergence ####
if (isTRUE(!lavaan::lavInspect(model.fit$model.fit, what = "converged"))) { stop("CFA model did not converge.", call. = FALSE) }
#—————————————————————————————————————— #
### Standard Error ####
if (isTRUE(se != "none")) { if (isTRUE(any(is.na(unlist(lavaan::lavInspect(model.fit$model.fit, what = "se")))))) { stop("Standard errors could not be computed.", call. = FALSE) } }
#—————————————————————————————————————— #
### Variance-Covariance Matrix of the Estimated Parameters ####
if (isTRUE(se != "none")) {
eigvals <- eigen(lavaan::lavInspect(model.fit$model.fit, what = "vcov"), symmetric = TRUE, only.values = TRUE)$values
# Model contains equality constraints
model.fit.par <- lavaan::parameterTable(model.fit$model.fit)$op == "=="
if (isTRUE(any(model.fit.par))) { eigvals <- rev(eigvals)[-seq_len(sum(model.fit.par))] }
if (isTRUE(min(eigvals) < .Machine$double.eps^(3L/4L))) {
warning("The variance-covariance matrix of the estimated parameters is not positive definite. This may be a symptom that the model is not identified.", call. = FALSE)
}
}
#—————————————————————————————————————— #
### Negative Variance of Observed Variables ####
#···················
#### Within Level
if (isTRUE(any(diag(lavaan::lavInspect(model.fit$model.fit, what = "theta")$within) < 0L))) {
warning("Some estimated variances of the observed variables at the Within level are negative.", call. = FALSE)
} else if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "theta")$within, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the residuals of the observed variables is not positive definite.", call. = FALSE)
}
#···················
#### Between Level
if (isTRUE(any(diag(lavaan::lavInspect(model.fit$model.fit, what = "theta")$.cluster) < 0L))) {
warning("Some estimated variances of the observed variables at the Between level are negative.", call. = FALSE)
} else if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "theta")$.cluster, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the residuals of the observed variables at the Between level is not positive definite.", call. = FALSE)
}
#—————————————————————————————————————— #
### Negative Variance of Latent Variables ####
#···················
#### Within Level
# Negative variance estimates
if (isTRUE(!is.null(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within))) {
if (isTRUE(any(diag(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within) < 0L))) {
warning("Some estimated variances of the latent variables at the Within level are negative.", call. = FALSE)
}
# Model-implied variance-covariance matrix of the latent variables
} else if (any(dim(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within) != 0L)) {
if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the latent variables at the Within level is not positive definite.", call. = FALSE)
}
}
#···················
#### Between Level
# Negative variance estimates
if (isTRUE(!is.null(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$cluster))) {
if (isTRUE(any(diag(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster) < 0L))) {
warning("Some estimated variances of the latent variables at the Between level are negative.", call. = FALSE)
}
# Model-implied variance-covariance matrix of the latent variables
} else if (any(dim(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster) != 0L)) {
if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the latent variables at the Between level is not positive definite.", call. = FALSE)
}
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Multilevel Reliability ####
#—————————————————————————————————————— #
### Within ####
# Within model parameter
param.w <- model.fit$result$param$within
# Within factor loading
load.w <- na.omit(param.w[param.w$param == "latent variable", "est"])
# Within residual variances
resid.w <- param.w[param.w$param == "residual variance", "est"]
# Fix negative residual variances at 0
resid.w <- ifelse(resid.w < 0L, 0L, resid.w)
# Within residual covariances
rescov.w <- na.omit(param.w[param.w$param == "residual covariance", "est"])
#—————————————————————————————————————— #
### Between ####
# Between model parameter
param.b <- model.fit$result$param$between
# Between factor loading
load.b <- na.omit(param.b[param.b$param == "latent variable", "est"])
# Between residual variances
resid.b <- param.b[param.b$param == "residual variance", "est"]
# Fix negative residual variances at 0
resid.b <- ifelse(resid.b < 0L, 0L, resid.b)
# Between factor variance
var.b <- param.b[param.b$param == "latent variance", "est"]
#—————————————————————————————————————— #
### Harmonic mean ####
hmean <- length(unique(model.fit$data$.cluster)) / sum(1L / table(x$.cluster) |> (\(p) p[p > 0L])())
#—————————————————————————————————————— #
### Level-Specific Omegas ####
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# Omega Within
omega.w <- sum(load.w)^2L / (sum(load.w)^2L + sum(resid.w) + 2L*sum(rescov.w))
#···················
#### Shared Cluster-Level Construct ####
}, shared = {
# Omega Between
omega.b <- sum(load.b)^2L / (sum(load.b)^2L + sum(resid.b) + ((sum(resid.w) + 2L*sum(rescov.w)) / hmean))
#···················
#### Configural Construct ####
}, config = {
# Omega Within
omega.w <- sum(load.w)^2L / (sum(load.w)^2L + sum(resid.w) + 2L*sum(rescov.w))
# Omega Between
omega.b <- (sum(load.b)^2L * var.b) / (sum(load.b)^2L * (1L / hmean + var.b) + sum(resid.b) + ((sum(resid.w) + 2L*sum(rescov.w)) / hmean))
# Overall Omega
omega.2l <- (sum(load.b)^2L * (1L + var.b)) / (sum(load.b)^2L * (1L + var.b) + sum(resid.b) + sum(resid.w) + 2L*sum(rescov.w))
})
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Confidence Interval for the Reliability ####
omega.2l.sim <- omega.b.sim <- omega.w.sim <- NULL
if (isTRUE(se != "none")) {
#—————————————————————————————————————— #
### Parameter Names ####
#···················
#### Within ####
# Within factor loading
load.w.par <- paste0("L", seq_len(length(load.w)))
# Within residual variances
resid.w.par <- apply(param.w[param.w$param == "residual variance", c("lhs", "op", "rhs")], 1L, paste, collapse = "")
# Within residual covariances
rescov.w.par <- apply(param.w[param.w$param == "residual covariance" & !is.na(param.w$est), c("lhs", "op", "rhs")], 1L, paste, collapse = "")
#···················
#### Between ####
# Between factor loading
load.b.par <- paste0("L", seq_len(length(load.b)))
# Between residual variances
resid.b.par <- paste0(apply(param.b[param.b$param == "residual variance", c("lhs", "op", "rhs")], 1L, paste, collapse = ""), ".l2")
# Between factor variance
var.b.par <- paste0(apply(param.b[param.b$param == "latent variance", c("lhs", "op", "rhs")], 1L, paste, collapse = ""), ".l2")
#···················
#### Within and Between ####
switch(const, within = {
parname <- c(load.w.par, resid.w.par, rescov.w.par)
}, shared = {
parname <- c(load.b.par, resid.b.par, resid.w.par, rescov.w.par)
}, config = {
parname <- c(load.w.par, load.b.par, var.b.par, resid.w.par, resid.b.par, rescov.w.par)
})
#—————————————————————————————————————— #
### Monte Carlo Confidence Interval ####
# Parameter estimates
fit.est <- na.omit(lavaan::coef(model.fit$model.fit)[parname])
# Variance-covariance matrix
fit.vcov <- lavaan::lavInspect(model.fit$model.fit, what = "vcov")[names(fit.est), names(fit.est)]
# Set seed
if (isTRUE(!is.null(seed))) { set.seed(seed) }
# Simulate from a multivariate normal distribution
simdata <- .internal.mvrnorm(nrep, mu = fit.est, Sigma = fit.vcov)
# Adapt parameter names
resid.w.par <- gsub("~~", "..", resid.w.par)
resid.b.par <- gsub("~~", "..", resid.b.par)
rescov.w.par <- gsub("~~", "..", rescov.w.par)
var.b.par <- gsub("~~", "..", var.b.par)
# Remove fixed Between residuals
if (isTRUE(!is.null(fix.resid))) {
resid.b.par <- resid.b.par[which(!resid.b.par %in% sapply(fix.resid, function(y) paste0(y, "..", y, ".l2")))]
# All Between residuals fixed at 0
if (isTRUE(length(resid.b.par) == 0L)) {
resid.b.par <- ".resid.b"
simdata <- data.frame(simdata, .resid.b = 0L)
}
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Multilevel Reliability Result Table ####
if (isTRUE(se != "none")) {
#—————————————————————————————————————— #
### Standard Errors and Confidence Intervals ####
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# No residual covariances at the Within level
if (isTRUE(length(rescov.w.par) == 0L)) {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), "))")))
} else {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")))")))
}
# Result table
omega <- data.frame(type = "omega.w", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.w,
low = quantile(omega.w.sim, probs = (1L - conf.level) / 2L),
upp = quantile(omega.w.sim, probs = 1L - (1L - conf.level) / 2L), row.names = NULL)
#···················
#### Shared-Cluster Construct ####
}, shared = {
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " + ", paste0(resid.b.par, collapse = " + "), " + (", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")) / hmean)", ")")))
# Result table
omega <- data.frame(type = "omega.b", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.b,
low = quantile(omega.b.sim, probs = (1L - conf.level) / 2L),
upp = quantile(omega.b.sim, probs = 1L - (1L - conf.level) / 2L), row.names = NULL)
#···················
#### Configural Construct ####
}, config = {
# No residual covariances at the Within level
if (isTRUE(length(rescov.w.par) == 0L)) {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), "))")))
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2 * ", var.b.par, ") / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " * (1 / ", hmean, " + ", var.b.par, ") + ", paste0(resid.b.par, collapse = " + "), " + ((", paste0(resid.w.par, collapse = " + "), ") / hmean)", "))"))))
# Overall Omega
eval(parse(text = paste0("omega.2l.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2L * (1L + var.b)) / ((", paste0(load.b.par, collapse = " + "), ")^2 * (1L + var.b) + ", paste0(resid.b.par, collapse = " + "), " + ", paste0(resid.w.par, collapse = " + "), "))"))))
} else {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), " + 2*( ", paste0(rescov.w.par, collapse = " + "), ")))")))
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2 * ", var.b.par, ") / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " * (1 / ", hmean, " + ", var.b.par, ") + ", paste0(resid.b.par, collapse = " + "), " + (", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")) / hmean)", ")"))))
# Overall Omega
eval(parse(text = paste0("omega.2l.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2L * (1L + var.b)) / ((", paste0(load.b.par, collapse = " + "), ")^2 * (1L + var.b) + ", paste0(resid.b.par, collapse = " + "), " + ", paste0(resid.w.par, collapse = " + "), " + 2L*(", paste0(rescov.w.par, collapse = " + "), ")))"))))
}
# Result table
omega <- data.frame(type = c("omega.w", "omega.b", "omega.2l"),
n.items = length(lavaan::lavNames(model.fit$model.fit)),
omega = c(omega.w, omega.b, omega.2l),
low = c(quantile(omega.w.sim, probs = (1L - conf.level) / 2L),
quantile(omega.b.sim, probs = (1L - conf.level) / 2L),
quantile(omega.2l.sim, probs = (1L - conf.level) / 2L)),
upp = c(quantile(omega.w.sim, probs = 1L - (1L - conf.level) / 2L),
quantile(omega.b.sim, probs = 1L - (1L - conf.level) / 2L),
quantile(omega.2l.sim, probs = 1L - (1L - conf.level) / 2L)), row.names = NULL)
})
#—————————————————————————————————————— #
### No Standard Errors and Confidence Intervals ####
} else {
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# Result table
omega <- data.frame(type = "omega.w", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.w, row.names = NULL)
#···················
#### Shared Cluster-Level Construct ####
}, shared = {
# Result table
omega <- data.frame(type = "omega.b", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.b)
#···················
#### Configural Construct ####
}, config = {
# Result table
omega <- data.frame(type = c("omega.w", "omega.b", "omega.2l"),
n.items = length(lavaan::lavNames(model.fit$model.fit)),
omega = c(omega.w, omega.b, omega.2l), row.names = NULL)
})
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Descriptive Statistics and Intraclass Correlation Coefficient, ICC(1) ####
itemstat <- NULL
if (isTRUE("item" %in% print)) {
switch(const,
within = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), wstd.ld = na.omit(param.w[param.w$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
}, shared = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), bstd.ld = na.omit(param.b[param.b$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
}, config = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), wstd.ld = na.omit(param.w[param.w$param == "latent variable", "stdyx"]), bstd.ld = na.omit(param.b[param.b$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
})
}
#_____________________________________________________________________________
#
# Return object --------------------------------------------------------------
object <- list(call = match.call(),
type = "multilevel.omega",
data = x,
args = list(rescov = rescov, const = const, fix.resid = fix.resid, se = se, optim.method = optim.method, missing = missing, nrep = nrep, seed = seed, conf.level = conf.level, print = print, digits = digits, r.digits = r.digits, as.na = as.na, write = write, append = append, check = check, output = output),
model = model.fit$model,
model.fit = model.fit$model.fit,
result = list(omega = omega, itemstat = itemstat))
class(object) <- "misty.object"
#_____________________________________________________________________________
#
# Write Results --------------------------------------------------------------
if (isTRUE(!is.null(write))) { .write.result(object = object, write = write, append = append) }
#_____________________________________________________________________________
#
# Output ---------------------------------------------------------------------
if (isTRUE(output)) { print(object, check = FALSE) }
return(invisible(object))
}
#_______________________________________________________________________________
Try the misty package in your browser
Any scripts or data that you put into this service are public.
misty documentation built on Aug. 2, 2026, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions multilevel.omega
Documented in multilevel.omega
#' Multilevel Composite Reliability
#'
#' This function computes point estimate and Monte Carlo confidence interval for
#' the multilevel composite reliability defined by Lai (2021) for a (1) within-cluster
#' construct, (2) shared cluster-level construct, and (3) individual and configural
#' cluster construct by calling the \code{cfa} function in the R package \pkg{lavaan}.
#' By default, the function prints level-specific multilevel composite reliability
#' indices for an individual and configural cluster construct with 95% Monte Carlo
#' confidence interval based on Huber-White standard errors.
#'
#' @param data a data frame. Multilevel confirmatory factor analysis
#' based on a measurement model with one factor at the Within
#' level and one factor at the Between level comprising all
#' variables in the data frame is conducted. Note that the
#' cluster variable specified in \code{cluster} is excluded
#' from \code{data} when specifying the argument \code{cluster}
#' using the variable name of the cluster variable.
#' @param ... an expression indicating the variable names in \code{data},
#' e.g., \code{multilevel.omega(dat, x1, x2, x3, cluster = "cluster)}.
#' Note that the operators \code{+}, \code{-},
#' \code{~}, \code{:}, \code{::}, and \code{!} can also be
#' used to select variables, see 'Details' in the
#' \code{\link{df.subset}} function.
#' @param cluster either a character string indicating the variable name of
#' the cluster variable in \code{data}, or a vector representing
#' the nested grouping structure (i.e., group or cluster variable).
#' @param rescov a character vector or a list of character vectors for specifying
#' residual covariances at the Within level, e.g. \code{rescov = c("x1", "x2")}
#' for specifying a residual covariance between indicators \code{x1}
#' and \code{x2} at the Within level or \code{rescov = list(c("x1", "x2"), c("x3", "x4"))}
#' for specifying residual covariances between indicators \code{x1}
#' and \code{x2}, and indicators \code{x3} and \code{x4} at
#' the Within level. Note that residual covariances at the
#' Between level cannot be specified using this function.
#' @param const a character string indicating the type of construct(s), i.e.,
#' \code{"within"} for within-cluster constructs, \code{"shared"}
#' for shared cluster-level constructs, and \code{"config"}
#' (default) for individual and configural constructs.
#' @param fix.resid a character vector for specifying residual variances to be
#' fixed at 0 at the Between level, e.g., \code{fix.resid = c("x1", "x3")}
#' to fix residual variances of indicators \code{x1} and \code{x2}
#' at the Between level at 0. Note that it is also possible
#' to specify \code{fix.resid = "all"} which fixes all residual
#' variances at the Between level at 0 in line with the strong
#' factorial measurement invariance assumption across cluster.
#' @param se a character string indicating the standard errors used for
#' computing Monte Carlo confidence intervals, i.e.,
#' \code{"none"} for no standard errors, \code{"standard"} for
#' conventional standard error based on inverting the expected
#' observed or first.order information matrix, and
#' \code{"robust.huber.white"} for the 'MLR' (aka pseudo ML,
#' Huber-White) approach. Note that \code{se = "none"} saves
#' computation time by not computing standard errors and
#' Monte Carlo confidence intervals.
#' @param optim.method a character string indicating the optimizer, i.e., \code{"nlminb"}
#' (default) for the unconstrained and bounds-constrained
#' quasi-Newton method optimizer and \code{"em"} for the
#' Expectation Maximization (EM) algorithm.
#' @param missing a character string indicating how to deal with missing data,
#' i.e., \code{"listwise"} for listwise deletion or \code{"fiml"}
#' (default) for full information maximum likelihood (FIML)
#' method.
#' @param nrep an integer value indicating the number of Monte Carlo
#' repetitions for computing confidence intervals.
#' @param seed a numeric value specifying the seed of the random number
#' generator for computing the Monte Carlo confidence interval.
#' @param conf.level a numeric value between 0 and 1 indicating the confidence
#' level of the interval.
#' @param print a character vector indicating which results to show, i.e.
#' \code{"all"} (default) for all results \code{"omega"} for
#' the composite reliability omega and \code{"item"} for item
#' statistics.
#' @param digits an integer value indicating the number of decimal places
#' to be used for displaying mean, standard deviation, minimum,
#' maximum, skewness, and kurtosis.
#' @param r.digits an integer value indicating the number of decimal places
#' to be used for displaying multilevel coefficient omega,
#' ICC(1), and standardized factor loadings.
#' @param as.na a numeric vector indicating user-defined missing values,
#' i.e. these values are converted to \code{NA} before conducting
#' the analysis. Note that \code{as.na()} function is only
#' applied to \code{data} but not to \code{cluster}.
#' @param write a character string naming a file for writing the output into
#' either a text file with file extension \code{".txt"} (e.g.,
#' \code{"Output.txt"}) or Excel file with file extension
#' \code{".xlsx"} (e.g., \code{"Output.xlsx"}). If the file
#' name does not contain any file extension, an Excel file will
#' be written.
#' @param append logical: if \code{TRUE} (default), output will be appended
#' to an existing text file with extension \code{.txt} specified
#' in \code{write}, if \code{FALSE} existing text file will be
#' overwritten.
#' @param check logical: if \code{TRUE} (default), argument specification,
#' convergence and model identification is checked.
#' @param output logical: if \code{TRUE} (default), output is shown.
#'
#' @details
#' Geldhof et al. (2014) introduced a multilevel confirmatory factor analysis
#' approach for estimating score reliability at each level. This approach uses
#' the estimated factor loadings, factor variances, and residual variances from both
#' the within-level and between-level model to compute the within-level reliability
#' \eqn{\tilde{\omega}^w} and the between-level reliability \eqn{\tilde{\omega}^b}
#' based on McDonald's \eqn{\omega}. However, Lai (2021) pointed out two limitations
#' of this approach. First, it does not account for different conceptualizations
#' of constructs in multilevel data, where constructs may have different meanings
#' (Stapleton et al., 2016). Second, the between-level reliability \eqn{\tilde{\omega}^b}
#' in Geldhof et al. (2014) is a measure of reliability of the latent cluster means,
#' which ignores the sampling error in the observed cluster means (see Lüdtke et
#' al., 2011). Consequently, \eqn{\tilde{\omega}^b} substantially overestimates
#' the true reliability of observed between-level composite scores. To address
#' these limitations, Lai (2021) proposed multilevel reliability coefficients
#' tailored to different types of constructs to estimate the reliability of the
#' observed scores at each level: (1) within-cluster construct, (2) shared
#' cluster-level construct, and (3) individual and configural constructs.
#'
#' \describe{
#' \item{\strong{Within-Cluster Construct}}{The latent variable representing the
#' within-cluster construct is only meaningful at the within level. For example,
#' students' popularity among their peers within the same classroom based on
#' sociometric ratings. In this case any variation across classrooms should be
#' irrelevant to the construct. Thus, a saturated between-level model is specified.
#' The reliability for a within-level composite, \eqn{Z^w_{ij}}, measuring
#' a within-cluster construct is
#'
#' \deqn{\omega^w = \frac{ (\sum^p_{k = 1}\lambda^w_k)^2\phi^w} {(\sum^p_{k = 1}\lambda^w_k)^2\phi^w + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' Note that the formula is the same as the formula for the individual construct
#' of an individual and configural construct, except that the factor loadings are
#' specific to the within level as the between-level model is saturated.}
#'
#' \item{\strong{Shared Cluster-Level Construct}}{The latent variable representing
#' the shared construct is only meaningful at the between level. For example,
#' safety of a neighborhood, effectiveness of a teacher, or organizational
#' climate are measured by multiple informants at the within level. In these cases
#' any variation within clusters should be irrelevant to the construct. Thus, a
#' saturated within-level model is specified. The reliability for a
#' between-level composite, \eqn{Z^b_{j}}, measuring a shared construct is
#'
#' \deqn{\omega^b = \frac{(\sum^p_{k = 1}\lambda^b_k)^2\phi^b} {(\sum^p_{k = 1}\lambda^b_k)^2\phi^b + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1}} + \mathbf{1}^{\prime}\mathbf{\sum}^w\mathbf{1}/\tilde{n}}
#'
#' where \eqn{\mathbf{\sum}^w} is the variance covariance matrix of the saturated
#' within-level model and \eqn{\tilde{n}} is the harmonic mean of cluster sizes
#' \eqn{\tilde{n} = \frac{1}{\sum^p_{k = 1}\frac{1}{n_j}}}.}
#'
#' \item{\strong{Individual and Configural Construct}}{The individual construct is
#' defined at the within level (e.g., individual student achievement), whereas the
#' configural construct consists of cluster averages of individual constructs (e.g.,
#' mean student achievement of a classroom). Note that the factor loadings need to
#' be constrained to be equal across levels, i.e., \eqn{\mathbf{\lambda}^w = \mathbf{\lambda}^b}
#' so that the latent variables at the within level, \eqn{\eta^w}, and between level,
#' \eqn{\eta^b}, are on the same metric.
#'
#' \describe{
#' \item{\emph{Individual Construct}}{The reliability for a within-level
#' composite, \eqn{Z^w_{ij}}, measuring an individual construct is
#'
#' \deqn{\omega^w = \frac{ (\sum^p_{k = 1}\lambda^k)^2\phi^w} {(\sum^p_{k = 1}\lambda^k)^2\phi^w + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' }
#'
#' \item{\emph{Configural Construct}}{The reliability for a between-level
#' composite, \eqn{Z^b_{j}}, measuring a configural construct is
#'
#' \deqn{\omega^b = \frac{(\sum^p_{k = 1}\lambda_k)^2\phi^b}{(\sum^p_{k = 1}\lambda_k)^2(\phi^b + \phi^w / \tilde{n}) + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1} + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}/\tilde{n}}}}
#'
#' \item{\emph{Composite Measuring an Individual and Configural Construct}}{
#' The reliability for an overall composite, \eqn{Z_{ij}}, measuring an individual
#' construct capturing the population variance of both the within-level and
#' the between-level components of the true score and the errors is
#'
#' \deqn{\omega^{2l} = \frac{(\sum^p_{k = 1}\lambda_k)^2(\phi^w + \phi^b)}{(\sum^p_{k = 1}\lambda_k)^2(\phi^w + \phi^b) + \mathbf{1}^{\prime}\mathbf{\Theta}^b\mathbf{1} + \mathbf{1}^{\prime}\mathbf{\Theta}^w\mathbf{1}}}
#'
#' Note that \eqn{\omega^{2l}} is simply the ratio of true variance to total
#' variance of the observed scores. In practice, \eqn{\omega^{2l}} is not
#' the most interesting coefficient because researchers are mainly interested
#' in distinguishing within-and between effects (Castro-Alvarez et al., 2026).}
#'
#' }
#' }
#' }
#' Note that the multilevel factor model assumes \eqn{p} items measuring one latent
#' construct at the within level (within-cluster construct), between level (shared
#' cluster-level construct), or within- and between level (individual and configural
#' construct). Local independence is assumed, \eqn{\Theta^w_j = diag[\theta^w_{11} ... \theta^w_{pp}]}
#' and \eqn{\Theta^b_j = diag[\theta^b_{11} ... \theta^b_{pp}]}, but composite reliability
#' can be computed when this assumption is violated by specifying the \code{rescov}
#' argument. The model assumes equal factor loadings across clusters,
#' \eqn{\mathbf{\lambda}^w_j = \mathbf{\lambda}^w}, and homogeneity of error
#' covariances across cluster, \eqn{\mathbf{\Theta}^w_j = \mathbf{\Theta}^w}.
#' In case of an individual and configural construct, the model also assumes cross-level
#' measurement invariance, i.e., \eqn{\mathbf{\lambda}^w = \mathbf{\lambda}^b}}.
#'
#' @author
#' Takuya Yanagida \email{takuya.yanagida@@univie.ac.at}
#'
#' @seealso
#' \code{\link{multilevel.alpha}}, \code{\link{item.omega}}, \code{\link{multilevel.cfa}},
#' \code{\link{multilevel.fit}}, \code{\link{multilevel.invar}}, \code{\link{multilevel.cor}},
#' \code{\link{multilevel.descript}}, \code{\link{write.result}}
#'
#' @references
#' Castro-Alvarez, S., Bringmann, L. F., Back, J., & Liu, S. (2026). The many
#' reliabilities of psychological dynamics: An overview of statistical approaches
#' to estimate the internal consistency reliability of intensive longitudinal data.
#' \emph{Psychological Methods, 31}(2), 281–296. https://doi.org/10.1037/met0000778
#'
#' Geldhof, G. J., Preacher, K. J., & Zyphur, M. J. (2014). Reliability
#' estimation in a multilevel confirmatory factor analysis framework. \emph{Psychological
#' Methods, 19}, 72-91. http://dx.doi.org/10.1037/a0032138
#'
#' Lai, M. H. C. (2021). Composite reliability of multilevel data: It’s about
#' observed scores and construct meanings. \emph{Psychological Methods, 26}(1),
#' 90–102. https://doi.org/10.1037/met0000287
#'
#' Lüdtke, O., Marsh, H. W., Robitzsch, A., & Trautwein, U. (2011). A 2 2
#' taxonomy of multilevel latent contextual models: Accuracy bias tradeoffs
#' in full and partial error correction models. \emph{Psychological Methods, 16},
#' 444-467. http://dx.doi.org/10.1037/a0024376
#'
#' Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling.
#' \emph{Journal of Statistical Software, 48}, 1-36. https://doi.org/10.18637/jss.v048.i02
#'
#' Stapleton, L. M., Yang, J. S., & Hancock, G. R. (2016). Construct meaning in
#' multilevel settings. \emph{Journal of Educational and Behavioral Statistics, 41}(5),
#' 481–520. https://doi.org/10.3102/1076998616646200
#'
#' Venables, W. N., Ripley, B. D. (2002).\emph{Modern Applied Statistics with S} (4th ed.).
#' Springer. https://www.stats.ox.ac.uk/pub/MASS4/.
#'
#' @return
#' Returns an object of class \code{misty.object}, which is a list with following
#' entries:
#' \item{\code{call}}{function call}
#' \item{\code{type}}{type of analysis}
#' \item{\code{data}}{data frame specified in \code{data} including the group variable
#' specified in \code{cluster}}
#' \item{\code{args}}{specification of function arguments}
#' \item{\code{model}}{specified model}
#' \item{\code{model.fit}}{fitted lavaan object (\code{mod.fit})}
#' \item{\code{check}}{results of the convergence and model identification check}
#' \item{\code{result}}{list with result tables, i.e., \code{omega} for the coefficient
#' omega including Monte Carlo confidence interval and
#' \code{itemstat} for descriptive statistics}
#'
#' @note
#' The function uses the functions \code{lavInspect}, \code{lavTech}, and \code{lavNames},
#' provided in the R package \pkg{lavaan} by Yves Rosseel (2012). The internal function
#' \code{.internal.mvrnorm} is a copy of the \code{mvrnorm} function in the package
#' \pkg{MASS} by Venables and Ripley (2002).
#'
#' @export
#'
#' @examples
#' \dontrun{
#' # Load data set "Demo.twolevel" in the lavaan package
#' data("Demo.twolevel", package = "lavaan")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Cluster Variable Specification
#'
#' # Example 1a: Specification using the argument '...'
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster")
#'
#' # Example 1b: Alternative specification with cluster variable 'cluster' in 'data'
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4", "cluster")], cluster = "cluster")
#'
#' # Example 1b: Alternative specification with cluster variable 'cluster' not in 'data'
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")], cluster = Demo.twolevel$cluster)
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Type of Construct
#'
#' # Example 2a: Within-Cluster Construct
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")],
#' cluster = Demo.twolevel$cluster, const = "within")
#'
#' # Example 2b: Shared Cluster-Level Construct
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "shared")
#'
#' # Example 2c: Configural Construct
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Residual (Co)Variances
#'
#' # Example 3a: Residual covariance between "y4" and "y5" at the Within level
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config",
#' rescov = c("y3", "y4"))
#'
#' # Example 3b: Residual variances of 'y1' at the Between level fixed at 0
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster", const = "config",
#' fix.resid = c("y1", "y2"))
#'
#' #————————————————————————————————————————————————————————————————————————————
# # Arguments 'se' and 'print'
#'
#' # Example 4a: No confidence intervals to speed up computation
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster", se = "none")
#'
#' # Example 4b: Omega and item statistics
#' multilevel.omega(Demo.twolevel, y1:y4, cluster = "cluster", print = "all")
#'
#' #————————————————————————————————————————————————————————————————————————————
#' # Write Results
#'
#' # Example 4a: Write results into a text file
#' multilevel.omega(Demo.twolevel[, c("y1", "y2", "y3", "y4")],
#' cluster = Demo.twolevel$cluster, write = "Multilevel_Omega.txt")
#'
#' # Example 4b: Write results into a Excel file
#' multilevel.omega(Demo.twolevel, y1, y2, y3, y4, cluster = "cluster",
#' write = "Multilevel_Omega.xlsx")
#' }
multilevel.omega <- function(data, ..., cluster, rescov = NULL,
const = c("within", "shared", "config"),
fix.resid = NULL, se = c("none", "standard", "robust.huber.white"),
optim.method = c("nlminb", "em"), missing = c("listwise", "fiml"),
nrep = 100000, seed = NULL, conf.level = 0.95,
print = c("all", "omega", "item"), digits = 2, r.digits = 3,
as.na = NULL, write = NULL, append = TRUE, check = TRUE,
output = TRUE) {
#_____________________________________________________________________________
#
# Initial Check --------------------------------------------------------------
# Check if input 'data' is missing or NULL
if (isTRUE(missing(data) || is.null(data))) { stop("Please specify a data frame for the argument 'data'", call. = FALSE) }
# Check input 'cluster'
if (isTRUE(missing(cluster) ||is.null(cluster))) { stop("Please specify a variable name or vector representing the grouping structure for the argument 'cluster'.", call. = FALSE) }
#_____________________________________________________________________________
#
# Data -----------------------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Using the Argument '...' ####
if (isTRUE(!missing(...))) {
# Extract data
x <- as.data.frame(data[, .var.names(data = data, ..., cluster = cluster), drop = FALSE])
# Extract cluster variable and convert tibble into data frame or vector
cluster <- data[, cluster] |> (\(y) if (isTRUE("tbl" %in% substr(class(y), 1L, 3L))) { unname(unlist(y)) } else { return(y) })()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Without Using the Argument '...' ####
} else {
# Data frame
x <- as.data.frame(data)
# Data and cluster
var.group <- .var.group(data = x, cluster = cluster)
# Data
if (isTRUE(!is.null(var.group$data))) { x <- var.group$data }
# Cluster variable
if (isTRUE(!is.null(var.group$cluster))) { cluster <- var.group$cluster }
}
#_____________________________________________________________________________
#
# Input Check ----------------------------------------------------------------
# Check inputs
.check.input(logical = c("append", "output"),
s.character = list(const = c("within", "shared", "config"), se = c("none", "standard", "robust.huber.white"), optim.method = c("nlminb", "em"), missing = c("listwise", "fiml")),
m.character = list(print = c("all", "omega", "item")),
args = c("conf.level", "digits", "r.digits", "nrep", "seed", "write2"), envir = environment(), input.check = check)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Additional Checks
if (isTRUE(check)) {
#—————————————————————————————————————— #
### Check Input 'rescov' ####
if (isTRUE(!is.null(rescov))) {
# Two variables for each residual covariance
if (isTRUE(is.list(rescov) && any(sapply(rescov, length) != 2L))) { stop("Please specify a list of character vectors for the argument 'rescov', where each element has two variable names", call. = FALSE)
} else { if (isTRUE(length(rescov) != 2L)) { stop("Please specify a character vector with two variable names for the argument 'rescov'", call. = FALSE) } }
# Variable in 'data'
(!unique(unlist(rescov)) %in% colnames(x)) |> (\(y) if (isTRUE(any(y))) { stop(paste0("Variables specified in the argument 'rescov' were not found in 'data': ", paste(unique(unlist(rescov))[y], collapse = ", ")), call. = FALSE) })()
}
#—————————————————————————————————————— #
### Check Input 'fix.resid' ####
(!unique(fix.resid) %in% colnames(x)) |> (\(y) if (isTRUE(any(y) && all(fix.resid != "all"))) { stop(paste0("Variables specified in the argument 'fix.resid' were not found in 'data': ", paste(fix.resid[y], collapse = ", ")), call. = FALSE) })()
}
#_____________________________________________________________________________
#
# Data and Arguments ---------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Manifest Variables ####
var <- colnames(x)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Data Frame with Cluster Variable ####
x <- data.frame(x, .cluster = cluster, row.names = NULL)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing Data on the Cluster Variable ####
if (isTRUE(any(is.na(x$.cluster)))) {
warning(paste0("Data contains missing values on the cluster variable, number of cases removed from the analysis: ", sum(is.na(x$.cluster))), call. = FALSE)
x <- x[!is.na(x$.cluster), ]
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Convert User-Missing Values into NA ####
if (isTRUE(!is.null(as.na))) { x[, var] <- .as.na(x[, var], na = as.na) }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Type of Construct(s) ####
if (isTRUE(all(c("within", "shared", "config") %in% const))) { const <- "config" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Residual Variances Fixed at 0 ####
if (isTRUE(fix.resid == "all")) { fix.resid <- var }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Standard Error ####
if (isTRUE(all(c("none", "standard", "robust.huber.white") %in% se))) { se <- "robust.huber.white" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Optimizer ####
if (isTRUE(all(c("nlminb", "em") %in% optim.method))) { optim.method <- "nlminb" }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing ####
# Complete data
if (isTRUE(all(!is.na(x[, var])))) {
missing <- "listwise"
# Data with missing values
} else {
if (isTRUE(all(c("listwise", "fiml") %in% missing))) {
missing <- "fiml"
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Missing Data on All Variable ####
(misty::na.prop(x[, var], append = FALSE) == 1L) |> (\(y) if (isTRUE(any(y) && missing == "fiml")) { warning(paste0("Data contains cases with missing values on all variables, number of cases removed from the analysis: ", sum(y)), call. = FALSE) })()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Print ####
if (isTRUE(all(c("all", "omega", "item") %in% print))) {
print <- "omega"
} else if (isTRUE(all(print == "all"))) {
print <- c("omega", "item")
}
#_____________________________________________________________________________
#
# Main Function --------------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Model Estimation ####
model.fit <- tryCatch(suppressWarnings(misty::multilevel.cfa(x, cluster = ".cluster", model = NULL, rescov = rescov,
model.w = NULL, model.b = NULL, rescov.w = NULL, rescov.b = NULL,
const = const, fix.resid = fix.resid, ident = "var", ls.fit = FALSE,
estimator = "ML", optim.method = optim.method, test = "none", se = se,
missing = missing, output = FALSE, check = FALSE)),
error = function(y) {
if (isTRUE(missing == "fiml")) {
stop("There was an estimation problem in lavaan, switching to missing = \"listwise\" might solve the problem.", call. = FALSE)
} else {
stop("There was an estimation problem in lavaan, model could not be estimated.", call. = FALSE)
}})
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Convergence and Model Identification Checks ####
if (isTRUE(check)) {
#—————————————————————————————————————— #
### Model Convergence ####
if (isTRUE(!lavaan::lavInspect(model.fit$model.fit, what = "converged"))) { stop("CFA model did not converge.", call. = FALSE) }
#—————————————————————————————————————— #
### Standard Error ####
if (isTRUE(se != "none")) { if (isTRUE(any(is.na(unlist(lavaan::lavInspect(model.fit$model.fit, what = "se")))))) { stop("Standard errors could not be computed.", call. = FALSE) } }
#—————————————————————————————————————— #
### Variance-Covariance Matrix of the Estimated Parameters ####
if (isTRUE(se != "none")) {
eigvals <- eigen(lavaan::lavInspect(model.fit$model.fit, what = "vcov"), symmetric = TRUE, only.values = TRUE)$values
# Model contains equality constraints
model.fit.par <- lavaan::parameterTable(model.fit$model.fit)$op == "=="
if (isTRUE(any(model.fit.par))) { eigvals <- rev(eigvals)[-seq_len(sum(model.fit.par))] }
if (isTRUE(min(eigvals) < .Machine$double.eps^(3L/4L))) {
warning("The variance-covariance matrix of the estimated parameters is not positive definite. This may be a symptom that the model is not identified.", call. = FALSE)
}
}
#—————————————————————————————————————— #
### Negative Variance of Observed Variables ####
#···················
#### Within Level
if (isTRUE(any(diag(lavaan::lavInspect(model.fit$model.fit, what = "theta")$within) < 0L))) {
warning("Some estimated variances of the observed variables at the Within level are negative.", call. = FALSE)
} else if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "theta")$within, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the residuals of the observed variables is not positive definite.", call. = FALSE)
}
#···················
#### Between Level
if (isTRUE(any(diag(lavaan::lavInspect(model.fit$model.fit, what = "theta")$.cluster) < 0L))) {
warning("Some estimated variances of the observed variables at the Between level are negative.", call. = FALSE)
} else if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "theta")$.cluster, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the residuals of the observed variables at the Between level is not positive definite.", call. = FALSE)
}
#—————————————————————————————————————— #
### Negative Variance of Latent Variables ####
#···················
#### Within Level
# Negative variance estimates
if (isTRUE(!is.null(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within))) {
if (isTRUE(any(diag(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within) < 0L))) {
warning("Some estimated variances of the latent variables at the Within level are negative.", call. = FALSE)
}
# Model-implied variance-covariance matrix of the latent variables
} else if (any(dim(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within) != 0L)) {
if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$within, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the latent variables at the Within level is not positive definite.", call. = FALSE)
}
}
#···················
#### Between Level
# Negative variance estimates
if (isTRUE(!is.null(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$cluster))) {
if (isTRUE(any(diag(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster) < 0L))) {
warning("Some estimated variances of the latent variables at the Between level are negative.", call. = FALSE)
}
# Model-implied variance-covariance matrix of the latent variables
} else if (any(dim(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster) != 0L)) {
if (isTRUE(any(eigen(lavaan::lavTech(model.fit$model.fit, what = "cov.lv")$.cluster, symmetric = TRUE, only.values = TRUE)$values < (-1L * .Machine$double.eps^(3/4))))) {
warning("The model-implied variance-covariance matrix of the latent variables at the Between level is not positive definite.", call. = FALSE)
}
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Multilevel Reliability ####
#—————————————————————————————————————— #
### Within ####
# Within model parameter
param.w <- model.fit$result$param$within
# Within factor loading
load.w <- na.omit(param.w[param.w$param == "latent variable", "est"])
# Within residual variances
resid.w <- param.w[param.w$param == "residual variance", "est"]
# Fix negative residual variances at 0
resid.w <- ifelse(resid.w < 0L, 0L, resid.w)
# Within residual covariances
rescov.w <- na.omit(param.w[param.w$param == "residual covariance", "est"])
#—————————————————————————————————————— #
### Between ####
# Between model parameter
param.b <- model.fit$result$param$between
# Between factor loading
load.b <- na.omit(param.b[param.b$param == "latent variable", "est"])
# Between residual variances
resid.b <- param.b[param.b$param == "residual variance", "est"]
# Fix negative residual variances at 0
resid.b <- ifelse(resid.b < 0L, 0L, resid.b)
# Between factor variance
var.b <- param.b[param.b$param == "latent variance", "est"]
#—————————————————————————————————————— #
### Harmonic mean ####
hmean <- length(unique(model.fit$data$.cluster)) / sum(1L / table(x$.cluster) |> (\(p) p[p > 0L])())
#—————————————————————————————————————— #
### Level-Specific Omegas ####
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# Omega Within
omega.w <- sum(load.w)^2L / (sum(load.w)^2L + sum(resid.w) + 2L*sum(rescov.w))
#···················
#### Shared Cluster-Level Construct ####
}, shared = {
# Omega Between
omega.b <- sum(load.b)^2L / (sum(load.b)^2L + sum(resid.b) + ((sum(resid.w) + 2L*sum(rescov.w)) / hmean))
#···················
#### Configural Construct ####
}, config = {
# Omega Within
omega.w <- sum(load.w)^2L / (sum(load.w)^2L + sum(resid.w) + 2L*sum(rescov.w))
# Omega Between
omega.b <- (sum(load.b)^2L * var.b) / (sum(load.b)^2L * (1L / hmean + var.b) + sum(resid.b) + ((sum(resid.w) + 2L*sum(rescov.w)) / hmean))
# Overall Omega
omega.2l <- (sum(load.b)^2L * (1L + var.b)) / (sum(load.b)^2L * (1L + var.b) + sum(resid.b) + sum(resid.w) + 2L*sum(rescov.w))
})
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Confidence Interval for the Reliability ####
omega.2l.sim <- omega.b.sim <- omega.w.sim <- NULL
if (isTRUE(se != "none")) {
#—————————————————————————————————————— #
### Parameter Names ####
#···················
#### Within ####
# Within factor loading
load.w.par <- paste0("L", seq_len(length(load.w)))
# Within residual variances
resid.w.par <- apply(param.w[param.w$param == "residual variance", c("lhs", "op", "rhs")], 1L, paste, collapse = "")
# Within residual covariances
rescov.w.par <- apply(param.w[param.w$param == "residual covariance" & !is.na(param.w$est), c("lhs", "op", "rhs")], 1L, paste, collapse = "")
#···················
#### Between ####
# Between factor loading
load.b.par <- paste0("L", seq_len(length(load.b)))
# Between residual variances
resid.b.par <- paste0(apply(param.b[param.b$param == "residual variance", c("lhs", "op", "rhs")], 1L, paste, collapse = ""), ".l2")
# Between factor variance
var.b.par <- paste0(apply(param.b[param.b$param == "latent variance", c("lhs", "op", "rhs")], 1L, paste, collapse = ""), ".l2")
#···················
#### Within and Between ####
switch(const, within = {
parname <- c(load.w.par, resid.w.par, rescov.w.par)
}, shared = {
parname <- c(load.b.par, resid.b.par, resid.w.par, rescov.w.par)
}, config = {
parname <- c(load.w.par, load.b.par, var.b.par, resid.w.par, resid.b.par, rescov.w.par)
})
#—————————————————————————————————————— #
### Monte Carlo Confidence Interval ####
# Parameter estimates
fit.est <- na.omit(lavaan::coef(model.fit$model.fit)[parname])
# Variance-covariance matrix
fit.vcov <- lavaan::lavInspect(model.fit$model.fit, what = "vcov")[names(fit.est), names(fit.est)]
# Set seed
if (isTRUE(!is.null(seed))) { set.seed(seed) }
# Simulate from a multivariate normal distribution
simdata <- .internal.mvrnorm(nrep, mu = fit.est, Sigma = fit.vcov)
# Adapt parameter names
resid.w.par <- gsub("~~", "..", resid.w.par)
resid.b.par <- gsub("~~", "..", resid.b.par)
rescov.w.par <- gsub("~~", "..", rescov.w.par)
var.b.par <- gsub("~~", "..", var.b.par)
# Remove fixed Between residuals
if (isTRUE(!is.null(fix.resid))) {
resid.b.par <- resid.b.par[which(!resid.b.par %in% sapply(fix.resid, function(y) paste0(y, "..", y, ".l2")))]
# All Between residuals fixed at 0
if (isTRUE(length(resid.b.par) == 0L)) {
resid.b.par <- ".resid.b"
simdata <- data.frame(simdata, .resid.b = 0L)
}
}
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Multilevel Reliability Result Table ####
if (isTRUE(se != "none")) {
#—————————————————————————————————————— #
### Standard Errors and Confidence Intervals ####
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# No residual covariances at the Within level
if (isTRUE(length(rescov.w.par) == 0L)) {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), "))")))
} else {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")))")))
}
# Result table
omega <- data.frame(type = "omega.w", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.w,
low = quantile(omega.w.sim, probs = (1L - conf.level) / 2L),
upp = quantile(omega.w.sim, probs = 1L - (1L - conf.level) / 2L), row.names = NULL)
#···················
#### Shared-Cluster Construct ####
}, shared = {
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " + ", paste0(resid.b.par, collapse = " + "), " + (", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")) / hmean)", ")")))
# Result table
omega <- data.frame(type = "omega.b", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.b,
low = quantile(omega.b.sim, probs = (1L - conf.level) / 2L),
upp = quantile(omega.b.sim, probs = 1L - (1L - conf.level) / 2L), row.names = NULL)
#···················
#### Configural Construct ####
}, config = {
# No residual covariances at the Within level
if (isTRUE(length(rescov.w.par) == 0L)) {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), "))")))
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2 * ", var.b.par, ") / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " * (1 / ", hmean, " + ", var.b.par, ") + ", paste0(resid.b.par, collapse = " + "), " + ((", paste0(resid.w.par, collapse = " + "), ") / hmean)", "))"))))
# Overall Omega
eval(parse(text = paste0("omega.2l.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2L * (1L + var.b)) / ((", paste0(load.b.par, collapse = " + "), ")^2 * (1L + var.b) + ", paste0(resid.b.par, collapse = " + "), " + ", paste0(resid.w.par, collapse = " + "), "))"))))
} else {
# Omega Within
eval(parse(text = paste0("omega.w.sim <- with(data.frame(simdata), ", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " / (", paste0("(", paste0(load.w.par, collapse = " + "), ")^2"), " + ", paste0(resid.w.par, collapse = " + "), " + 2*( ", paste0(rescov.w.par, collapse = " + "), ")))")))
# Omega Between
eval(parse(text = paste0("omega.b.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2 * ", var.b.par, ") / (", paste0("(", paste0(load.b.par, collapse = " + "), ")^2"), " * (1 / ", hmean, " + ", var.b.par, ") + ", paste0(resid.b.par, collapse = " + "), " + (", paste0(resid.w.par, collapse = " + "), " + 2*(", paste0(rescov.w.par, collapse = " + "), ")) / hmean)", ")"))))
# Overall Omega
eval(parse(text = paste0("omega.2l.sim <- with(data.frame(simdata), ", paste0("((", paste0(load.b.par, collapse = " + "), ")^2L * (1L + var.b)) / ((", paste0(load.b.par, collapse = " + "), ")^2 * (1L + var.b) + ", paste0(resid.b.par, collapse = " + "), " + ", paste0(resid.w.par, collapse = " + "), " + 2L*(", paste0(rescov.w.par, collapse = " + "), ")))"))))
}
# Result table
omega <- data.frame(type = c("omega.w", "omega.b", "omega.2l"),
n.items = length(lavaan::lavNames(model.fit$model.fit)),
omega = c(omega.w, omega.b, omega.2l),
low = c(quantile(omega.w.sim, probs = (1L - conf.level) / 2L),
quantile(omega.b.sim, probs = (1L - conf.level) / 2L),
quantile(omega.2l.sim, probs = (1L - conf.level) / 2L)),
upp = c(quantile(omega.w.sim, probs = 1L - (1L - conf.level) / 2L),
quantile(omega.b.sim, probs = 1L - (1L - conf.level) / 2L),
quantile(omega.2l.sim, probs = 1L - (1L - conf.level) / 2L)), row.names = NULL)
})
#—————————————————————————————————————— #
### No Standard Errors and Confidence Intervals ####
} else {
switch(const,
#···················
#### Within-Cluster Construct ####
within = {
# Result table
omega <- data.frame(type = "omega.w", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.w, row.names = NULL)
#···················
#### Shared Cluster-Level Construct ####
}, shared = {
# Result table
omega <- data.frame(type = "omega.b", n.items = length(lavaan::lavNames(model.fit$model.fit)), omega = omega.b)
#···················
#### Configural Construct ####
}, config = {
# Result table
omega <- data.frame(type = c("omega.w", "omega.b", "omega.2l"),
n.items = length(lavaan::lavNames(model.fit$model.fit)),
omega = c(omega.w, omega.b, omega.2l), row.names = NULL)
})
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Descriptive Statistics and Intraclass Correlation Coefficient, ICC(1) ####
itemstat <- NULL
if (isTRUE("item" %in% print)) {
switch(const,
within = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), wstd.ld = na.omit(param.w[param.w$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
}, shared = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), bstd.ld = na.omit(param.b[param.b$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
}, config = {
itemstat <- data.frame(misty::df.rename(model.fit$result$descript, from = "variable", to = "item"), wstd.ld = na.omit(param.w[param.w$param == "latent variable", "stdyx"]), bstd.ld = na.omit(param.b[param.b$param == "latent variable", "stdyx"]), fix.empty.names = FALSE)
})
}
#_____________________________________________________________________________
#
# Return object --------------------------------------------------------------
object <- list(call = match.call(),
type = "multilevel.omega",
data = x,
args = list(rescov = rescov, const = const, fix.resid = fix.resid, se = se, optim.method = optim.method, missing = missing, nrep = nrep, seed = seed, conf.level = conf.level, print = print, digits = digits, r.digits = r.digits, as.na = as.na, write = write, append = append, check = check, output = output),
model = model.fit$model,
model.fit = model.fit$model.fit,
result = list(omega = omega, itemstat = itemstat))
class(object) <- "misty.object"
#_____________________________________________________________________________
#
# Write Results --------------------------------------------------------------
if (isTRUE(!is.null(write))) { .write.result(object = object, write = write, append = append) }
#_____________________________________________________________________________
#
# Output ---------------------------------------------------------------------
if (isTRUE(output)) { print(object, check = FALSE) }
return(invisible(object))
}
#_______________________________________________________________________________
Try the misty package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.