Show factor correlations for oblique rotations

NAMESPACE

@@ -149,6 +149,8 @@ S3method(equivalence_test,rma)
 S3method(equivalence_test,wbm)
 S3method(equivalence_test,zeroinfl)
 S3method(factor_analysis,data.frame)
+S3method(factor_scores,fa)
+S3method(factor_scores,parameters_efa)
 S3method(format,compare_parameters)
 S3method(format,equivalence_test_lm)
 S3method(format,p_calibrate)
@@ -939,6 +941,7 @@ export(dominance_analysis)
 export(efa_to_cfa)
 export(equivalence_test)
 export(factor_analysis)
+export(factor_scores)
 export(format_df_adjust)
 export(format_order)
 export(format_p_adjust)

NEWS.md

@@ -9,6 +9,9 @@
 
 ## Changes
 
+* New function `factor_scores()` to extract factor scores from EFA (`psych::fa()`
+  or `factor_analysis()`).
+
 * Added and/or improved print-methods for all functions around PCA and FA.
 
 * Improved efficiency in `model_parameters()` for models from packages *brms*

R/display.R

@@ -172,7 +172,6 @@
       fun_args,
       list(
         column_labels = column_labels,
-        align = align,
         font_size = font_size,
         line_padding = line_padding
       )
@@ -237,7 +236,7 @@
 #' @inheritParams model_parameters.principal
 #' @rdname display.parameters_model
 #' @export
 display.parameters_efa <- function(object, format = "markdown", digits = 2, sort = FALSE, threshold = NULL, labels = NULL, ...) {
   format <- insight::validate_argument(format, c("markdown", "html", "md", "tt"))
 
   fun_args <- list(

R/factor_analysis.R

@@ -40,3 +40,8 @@ factor_analysis.data.frame <- function(x,
   attr(out, "dataset") <- x
   out
 }
+
+
+.is_oblique_rotation <- function(rotation) {
+  !is.null(rotation) && tolower(rotation) %in% c("promax", "oblimin", "simplimax", "bentlerQ", "geominQ", "biquartimin", "cluster") # nolint
+}

R/factor_scores.R

@@ -0,0 +1,34 @@
+#' Extract factor scores from Factor Analysis (EFA)
+#'
+#' `factor_scores()` extracts the factor scores from objects returned by
+#' [`psych::fa()`] or [`factor_analysis()`].
+#'
+#' @param x An object returned by [`psych::fa()`] or [`factor_analysis()`].
+#' @param ... Currently unused.
+#'
+#' @return A data frame with the factor scores. It simply extracts the `$scores`
+#' element from the object and converts it into a data frame.
+#'
+#' @examplesIf insight::check_if_installed("psych", quietly = TRUE)
+#' data(mtcars)
+#' out <- factor_analysis(mtcars[, 1:7], n = 2)
+#' head(factor_scores(out))
+#'
+#' @export
+factor_scores <- function(x, ...) {
+  UseMethod("factor_scores")
+}
+
+#' @export
+factor_scores.fa <- function(x, ...) {
+  as.data.frame(x$scores)
+}
+
+#' @export
+factor_scores.parameters_efa <- function(x, ...) {
+  model <- attributes(x)$model
+  if (is.null(model)) {
+    insight::format_error("The `model` attribute is missing from the input object.")
+  }
+  as.data.frame(model$scores)
+}

R/print_html.R

@@ -281,6 +281,18 @@ print_html.parameters_efa_summary <- function(x, digits = 3, ...) {
   } else if ("Component" %in% names(x)) {
     names(x) <- c("Component", "Eigenvalues", "Variance Explained", "Variance Explained (Cumulative)", "Variance Explained (Proportion)") # nolint
   }
+
+  # we may have factor correlations
+  fc <- attributes(x)$factor_correlations
+
+  # if we have factor correlations, we need to add them to the table
+  if (!is.null(fc)) {
+    fc$Component <- "Factor Correlations"
+    x$Component <- "Explained Variance"
+    colnames(fc)[1] <- colnames(x)[1]
+    x <- .safe(rbind(x, fc), x)
+  }
+
   insight::export_table(x, digits = digits, format = "html", caption = table_caption, align = "firstleft")
 }
 

R/print_md.R

@@ -254,6 +254,19 @@ print_md.parameters_efa_summary <- function(x, digits = 3, ...) {
   } else if ("Component" %in% names(x)) {
     names(x) <- c("Component", "Eigenvalues", "Variance Explained", "Variance Explained (Cumulative)", "Variance Explained (Proportion)") # nolint
   }
+
+  # we may have factor correlations
+  fc <- attributes(x)$factor_correlations
+
+  # if we have factor correlations, we need to add them to the table
+  if (!is.null(fc)) {
+    x <- list(x, fc)
+    table_caption <- list(
+      table_caption,
+      "Factor Correlations"
+    )
+  }
+
   insight::export_table(x, digits = digits, format = "markdown", caption = table_caption, align = "firstleft")
 }
 

R/utils_pca_efa.R

@@ -17,25 +17,25 @@
 #' score for each component from the PCA, which is on the same scale as the
 #' original, single items that were used to compute the PCA.
 #'
-#' @examples
-#' if (require("psych")) {
-#'   pca <- principal_components(mtcars[, 1:7], n = 2, rotation = "varimax")
+#' @return A data frame with subscales, which are average sum scores for all
+#'   items from each component.
 #'
-#'   # PCA extracted two components
-#'   pca
+#' @examplesIf insight::check_if_installed("psych", quietly = TRUE)
+#' pca <- principal_components(mtcars[, 1:7], n = 2, rotation = "varimax")
 #'
-#'   # assignment of items to each component
-#'   closest_component(pca)
+#' # PCA extracted two components
+#' pca
 #'
-#'   # now we want to have sum scores for each component
-#'   get_scores(pca)
+#' # assignment of items to each component
+#' closest_component(pca)
+#'
+#' # now we want to have sum scores for each component
+#' get_scores(pca)
+#'
+#' # compare to manually computed sum score for 2nd component, which
+#' # consists of items "hp" and "qsec"
+#' (mtcars$hp + mtcars$qsec) / 2
 #'
-#'   # compare to manually computed sum score for 2nd component, which
-#'   # consists of items "hp" and "qsec"
-#'   (mtcars$hp + mtcars$qsec) / 2
-#' }
-#' @return A data frame with subscales, which are average sum scores for all
-#'   items from each component.
 #' @export
 get_scores <- function(x, n_items = NULL) {
   subscales <- closest_component(x)
@@ -99,12 +99,21 @@ summary.parameters_efa <- function(object, ...) {
     colnames(x)
   )
 
-
   x <- as.data.frame(t(x[, cols]))
   x <- cbind(data.frame(Parameter = row.names(x), stringsAsFactors = FALSE), x)
   names(x) <- c("Parameter", attributes(object)$summary$Component)
   row.names(x) <- NULL
 
+  if (.is_oblique_rotation(attributes(object)$rotation)) {
+    factor_correlations <- attributes(object)$model$Phi
+    if (!is.null(factor_correlations)) {
+      attr(x, "factor_correlations") <- datawizard::rownames_as_column(
+        as.data.frame(factor_correlations),
+        var = "Factor"
+      )
+    }
+  }
+
   if (inherits(object, "parameters_efa")) {
     class(x) <- c("parameters_efa_summary", class(object))
   } else {
@@ -206,6 +215,9 @@ predict.parameters_pca <- predict.parameters_efa
 
 #' @export
 print.parameters_efa_summary <- function(x, digits = 3, ...) {
+  # we may have factor correlations
+  fc <- attributes(x)$factor_correlations
+
   if ("Parameter" %in% names(x)) {
     x$Parameter <- c(
       "Eigenvalues", "Variance Explained", "Variance Explained (Cumulative)",
@@ -225,6 +237,19 @@ print.parameters_efa_summary <- function(x, digits = 3, ...) {
     format = "text",
     ...
   ))
+
+  if (!is.null(fc)) {
+    cat("\n")
+    cat(insight::export_table(
+      fc,
+      digits = digits,
+      caption = c("# Factor Correlations", "blue"),
+      format = "text",
+      ...
+    ))
+  }
+
+
   invisible(x)
 }
 

pkgdown/_pkgdown.yml

@@ -85,6 +85,7 @@ reference:
     contents:
       - convert_efa_to_cfa
       - factor_analysis
+      - factor_scores
       - get_scores
       - n_factors
       - principal_components

tests/testthat/_snaps/factor_analysis.md

@@ -0,0 +1,44 @@
+# factor_analysis
+
+    Code
+      print(summary(out))
+    Output
+      # (Explained) Variance of Components
+
+      Parameter                       |   MR1 |   MR2
+      -----------------------------------------------
+      Eigenvalues                     | 4.947 | 1.062
+      Variance Explained              | 0.638 | 0.220
+      Variance Explained (Cumulative) | 0.638 | 0.858
+      Variance Explained (Proportion) | 0.744 | 0.256
+
+      # Factor Correlations
+
+      Factor |    MR1 |    MR2
+      ------------------------
+      MR1    |  1.000 | -0.366
+      MR2    | -0.366 |  1.000
+
+---
+
+    Code
+      print_md(summary(out))
+    Output
+
+
+      Table: (Explained) Variance of Components
+
+      |Parameter                       |   MR1 |   MR2 |
+      |:-------------------------------|:-----:|:-----:|
+      |Eigenvalues                     | 4.947 | 1.062 |
+      |Variance Explained              | 0.638 | 0.220 |
+      |Variance Explained (Cumulative) | 0.638 | 0.858 |
+      |Variance Explained (Proportion) | 0.744 | 0.256 |
+
+      Table: Factor Correlations
+
+      |Factor |    MR1 |    MR2 |
+      |:------|:------:|:------:|
+      |MR1    |  1.000 | -0.366 |
+      |MR2    | -0.366 |  1.000 |
+

tests/testthat/test-factor_analysis.R

@@ -1,4 +1,4 @@
-test_that("n_factors, default", {
+test_that("factor_analysis", {
   skip_on_cran()
   skip_if_not_installed("GPArotation")
   skip_if_not_installed("psych")
@@ -9,7 +9,7 @@ test_that("n_factors, default", {
   raq_items <- as.data.frame(discovr::raq)
   raq_items$id <- NULL
 
-  out <- parameters::factor_analysis(
+  out <- factor_analysis(
     raq_items,
     n = 4,
     scores = "tenBerge",
@@ -34,4 +34,19 @@ test_that("n_factors, default", {
     tolerance = 1e-3,
     ignore_attr = TRUE
   )
+
+  # include factor correlations
+  out <- factor_analysis(
+    mtcars[, 1:7],
+    n = 2,
+    rotation = "oblimin",
+    threshold = "max",
+    sort = TRUE
+  )
+  expect_snapshot(print(summary(out)))
+  expect_snapshot(print_md(summary(out)))
+
+  # check factor scores
+  fc <- factor_scores(out)
+  expect_identical(dim(fc), c(32L, 2L))
 })