draft

Author: strengejacke

NAMESPACE

@@ -138,6 +138,7 @@ S3method(logLik,plm)
 S3method(logLik,svycoxph)
 S3method(mcdonalds_omega,data.frame)
 S3method(mcdonalds_omega,matrix)
+S3method(mcdonalds_omega,parameters_pca)
 S3method(model_performance,Arima)
 S3method(model_performance,BFBayesFactor)
 S3method(model_performance,DirichletRegModel)
@@ -293,6 +294,7 @@ S3method(print,icc_decomposed)
 S3method(print,item_difficulty)
 S3method(print,item_discrimination)
 S3method(print,looic)
+S3method(print,mcdonalds_omega)
 S3method(print,performance_accuracy)
 S3method(print,performance_cv)
 S3method(print,performance_hosmer)

R/cronbachs_alpha.R

@@ -4,7 +4,8 @@
 #' @description Compute various measures of internal consistencies
 #'    for tests or item-scales of questionnaires.
 #'
-#' @param x A matrix or a data frame.
+#' @param x A matrix or a data frame, or an object returned by
+#' `[parameters::principal_components()]`.
 #' @param ... Currently not used.
 #'
 #' @return The Cronbach's Alpha value for `x`.
@@ -50,14 +51,12 @@ cronbachs_alpha.data.frame <- function(x, verbose = TRUE, ...) {
 }
 
 
-
 #' @export
 cronbachs_alpha.matrix <- function(x, verbose = TRUE, ...) {
   cronbachs_alpha(as.data.frame(x), verbose = verbose, ...)
 }
 
 
-
 #' @export
 cronbachs_alpha.parameters_pca <- function(x, verbose = TRUE, ...) {
   # fetch data used for the PCA

R/mcdonalds_omega.r

@@ -4,8 +4,9 @@
 #' @description Compute various measures of internal consistencies
 #'    for tests or item-scales of questionnaires.
 #'
-#' @param x A matrix or a data frame.
-#' @param ... Currently not used.
+#' @param ci Confidence interval for the reliability estimate. If `NULL`,
+#' no confidence interval is computed.
+#' @inheritParams cronbachs_alpha
 #'
 #' @return The McDonald's Omega value for `x`.
 #'
@@ -33,11 +34,24 @@ mcdonalds_omega <- function(x, ...) {
 
 #' @export
 mcdonalds_omega.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
-  varnames <- colnames(x)
-  n_params <- length(varnames)
-  name_loadings <- paste0("a", 1:n_params)
-  name_error <- paste0("b", 1:n_params)
+  # remove missings
+  .data <- stats::na.omit(x)
 
+  # we need at least two columns for Cronach's Alpha
+  if (is.null(ncol(.data)) || ncol(.data) < 2) {
+    if (verbose) {
+      insight::format_warning("Too few columns in `x` to compute McDonald's Omega.")
+    }
+    return(NULL)
+  }
+
+  # prepare names and formulas for lavaan
+  varnames <- colnames(.data)
+  name_loadings <- paste0("a", seq_len(ncol(.data)))
+  name_error <- paste0("b", seq_len(ncol(.data)))
+
+  # we need this specific formulation for lavaan to get the omega reliability estimate
+  # see code in MBESS
   model <- paste0("f1 =~ NA*", varnames[1], " + ")
   formula_loadings <- paste(paste0(name_loadings, "*", varnames), collapse = " + ")
   formula_factors <- "f1 ~~ 1*f1\n"
@@ -60,34 +74,102 @@ mcdonalds_omega.data.frame <- function(x, ci = 0.95, verbose = TRUE, ...) {
 
   insight::check_if_installed("lavaan")
 
-  fit <- lavaan::cfa(model, data = x, missing = "ml", estimator = "mlr", se = "default")
+  # fit CFA to get reliability estimate
+  fit <- lavaan::cfa(model, data = .data, missing = "ml", estimator = "mlr", se = "default")
   out <- lavaan::parameterEstimates(fit)
 
+  # extract omega and related standard error
   estimate <- as.vector(out$est[out$label == "relia"])
   se <- as.vector(out$se[out$label == "relia"])
 
+  # if user requested CI, return data frame with omega and CI
   if (!is.null(ci) && !is.na(ci)) {
     crit <- stats::qnorm((1 + ci) / 2)
 
     logest <- log(estimate / (1 - estimate))
     logse <- se / (estimate * (1 - estimate))
     loglower <- logest - crit * logse
     logupper <- logest + crit * logse
+
     if (logupper < loglower) {
       temp <- loglower
       loglower <- logupper
       loguppper <- temp
     }
-    ci_low <- 1 / (1 + exp(-loglower))
-    ci_high <- 1 / (1 + exp(-logupper))
+
+    omega <- data.frame(
+      Omega = estimate,
+      CI_low = 1 / (1 + exp(-loglower)),
+      CI_high = 1 / (1 + exp(-logupper)),
+      stringsAsFactors = FALSE
+    )
+    class(omega) <- c("mcdonalds_omega", "data.frame")
   } else {
-    ci_low <- NA
-    ci_high <- NA
+    omega <- estimate
   }
+
+  omega
 }
 
 
 #' @export
-mcdonalds_omega.matrix <- function(x, verbose = TRUE, ...) {
-  mcdonalds_omega(as.data.frame(x), verbose = verbose, ...)
+mcdonalds_omega.matrix <- function(x, ci = 0.95, verbose = TRUE, ...) {
+  mcdonalds_omega(as.data.frame(x), ci = ci, verbose = verbose, ...)
+}
+
+
+#' @export
+mcdonalds_omega.parameters_pca <- function(x, verbose = TRUE, ...) {
+  # fetch data used for the PCA
+  pca_data <- attributes(x)$dataset
+
+  # if NULL, can we get from environment?
+  if (is.null(pca_data)) {
+    pca_data <- attr(x, "data")
+    if (is.null(pca_data)) {
+      if (verbose) {
+        insight::format_warning("Could not find data frame that was used for the PCA.")
+      }
+      return(NULL)
+    }
+    pca_data <- get(pca_data, envir = parent.frame())
+  }
+
+  # get assignment of columns to extracted components, based on the max loading
+  factor_assignment <- attributes(x)$closest_component
+
+  # sort and get unique IDs so we only get data from relevant columns
+  unique_factors <- sort(unique(factor_assignment))
+
+  # apply cronbach's alpha for each component,
+  # only for variables with max loading
+  omegas <- sapply(unique_factors, function(i) {
+    mcdonalds_omega(
+      pca_data[, as.vector(x$Variable[factor_assignment == i]), drop = FALSE],
+      ci = NULL,
+      verbose = verbose,
+      ...
+    )
+  })
+
+  names(omegas) <- paste0("PC", unique_factors)
+  unlist(omegas)
+}
+
+
+# methods ---------------------------------------------------------------------
+
+#' @export
+print.mcdonalds_omega <- function(x, digits = 3, ...) {
+  # print regular R2
+  out <- sprintf(
+    "Omega: %.*f %s",
+    digits,
+    x$omega,
+    insight::format_ci(ci_low, ci_high, digits = digits, ci = NULL)
+  )
+
+  cat(out)
+  cat("\n")
+  invisible(x)
 }