diff --git a/NAMESPACE b/NAMESPACE
index 66529502a..ba98b0bf8 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index 3c0c83814..7ec62b0c6 100644
--- a/NEWS.md
+++ b/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*
diff --git a/R/display.R b/R/display.R
index 5fb158814..c35c80fd1 100644
--- a/R/display.R
+++ b/R/display.R
@@ -172,7 +172,6 @@ display.compare_parameters <- function(object,
       fun_args,
       list(
         column_labels = column_labels,
-        align = align,
         font_size = font_size,
         line_padding = line_padding
       )
diff --git a/R/factor_analysis.R b/R/factor_analysis.R
index ab4107061..423285296 100644
--- a/R/factor_analysis.R
+++ b/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
+}
diff --git a/R/factor_scores.R b/R/factor_scores.R
new file mode 100644
index 000000000..8a1648a42
--- /dev/null
+++ b/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)
+}
diff --git a/R/print_html.R b/R/print_html.R
index bef10671e..dfa8cfcc6 100644
--- a/R/print_html.R
+++ b/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")
 }
 
diff --git a/R/print_md.R b/R/print_md.R
index b7e6c9ee0..959a495c7 100644
--- a/R/print_md.R
+++ b/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")
 }
 
diff --git a/R/utils_pca_efa.R b/R/utils_pca_efa.R
index 24dd5305c..46c9f2f2c 100644
--- a/R/utils_pca_efa.R
+++ b/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)
 }
 
diff --git a/man/factor_scores.Rd b/man/factor_scores.Rd
new file mode 100644
index 000000000..8f8b2a5cf
--- /dev/null
+++ b/man/factor_scores.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/factor_scores.R
+\name{factor_scores}
+\alias{factor_scores}
+\title{Extract factor scores from Factor Analysis (EFA)}
+\usage{
+factor_scores(x, ...)
+}
+\arguments{
+\item{x}{An object returned by \code{\link[psych:fa]{psych::fa()}} or \code{\link[=factor_analysis]{factor_analysis()}}.}
+
+\item{...}{Currently unused.}
+}
+\value{
+A data frame with the factor scores. It simply extracts the \verb{$scores}
+element from the object and converts it into a data frame.
+}
+\description{
+\code{factor_scores()} extracts the factor scores from objects returned by
+\code{\link[psych:fa]{psych::fa()}} or \code{\link[=factor_analysis]{factor_analysis()}}.
+}
+\examples{
+\dontshow{if (insight::check_if_installed("psych", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+data(mtcars)
+out <- factor_analysis(mtcars[, 1:7], n = 2)
+head(factor_scores(out))
+\dontshow{\}) # examplesIf}
+}
diff --git a/man/get_scores.Rd b/man/get_scores.Rd
index 313e4519c..bc963728a 100644
--- a/man/get_scores.Rd
+++ b/man/get_scores.Rd
@@ -31,20 +31,20 @@ 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")
+\dontshow{if (insight::check_if_installed("psych", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+pca <- principal_components(mtcars[, 1:7], n = 2, rotation = "varimax")
 
-  # PCA extracted two components
-  pca
+# PCA extracted two components
+pca
 
-  # assignment of items to each component
-  closest_component(pca)
+# assignment of items to each component
+closest_component(pca)
 
-  # now we want to have sum scores for each component
-  get_scores(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
+\dontshow{\}) # examplesIf}
 }
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index d8dd28ca3..0d2752b65 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -85,6 +85,7 @@ reference:
     contents:
       - convert_efa_to_cfa
       - factor_analysis
+      - factor_scores
       - get_scores
       - n_factors
       - principal_components
diff --git a/tests/testthat/_snaps/factor_analysis.md b/tests/testthat/_snaps/factor_analysis.md
new file mode 100644
index 000000000..694d35d64
--- /dev/null
+++ b/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 |
+
diff --git a/tests/testthat/test-factor_analysis.R b/tests/testthat/test-factor_analysis.R
index 12e5f7d00..774dc5e17 100644
--- a/tests/testthat/test-factor_analysis.R
+++ b/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))
 })