jasp-stats · juliuspfadt · Jul 5, 2026 · Jul 5, 2026 · Jul 5, 2026 · Jul 5, 2026
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,10 +1,11 @@
 import(jaspBase)
 export(latentClassAnalysis)
 export(latentClassAnalysisInternal)
+export(numberOfFactors)
+export(numberOfFactorsInternal)
 export(principalComponentAnalysisInternal)
 export(exploratoryFactorAnalysisInternal)
 export(confirmatoryFactorAnalysisInternal)
 export(principalComponentAnalysis)
 export(exploratoryFactorAnalysis)
 export(confirmatoryFactorAnalysis)
-
diff --git a/NEWS.md b/NEWS.md
@@ -17,6 +17,13 @@
 
 ## Added
 * Latent Class Analysis ([PR #346](https://github.com/jasp-stats/jaspFactor/pull/346)).
+* New **Number of Factors/Components** analysis that consolidates all factor/component retention methods (parallel analysis, eigenvalue criterion, scree plot) into a dedicated analysis, usable for both EFA and PCA decisions.
+
+## Changed
+* EFA, PCA: eigenvalue- and parallel-analysis-based factor count options removed from both analyses; the number of factors/components is now set manually, with guidance to use the new Number of Factors/Components analysis.
+
+## Fixed
+* Number of Factors: eigenvalue criterion now uses principal component eigenvalues when PC-based parallel analysis is selected, and factor eigenvalues when FA-based parallel analysis is selected. Previously PC eigenvalues were always used, causing factors with factor eigenvalues below the threshold to be incorrectly retained in EFA ([jasp-issues#3970](https://github.com/jasp-stats/jasp-issues/issues/3970)).
 
 ## Fixed
 * EFA: polychoric/tetrachoric correlation matrix errors are now caught and shown as a user-friendly message, including which variables have missing response categories.

diff --git a/R/exploratoryfactoranalysis.R b/R/exploratoryfactoranalysis.R
@@ -28,8 +28,6 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
   if (ready)
     .pcaCheckErrors(dataset, options, method = "efa")
 
-  options(mc.cores = 1) # prevent the fa.parallel function using mutlitple cores by default
-
   modelContainer <- .efaModelContainer(jaspResults)
 
   # output functions
@@ -44,8 +42,6 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
   .efaCorrTable(            modelContainer, dataset, options, ready)
   .efaAdditionalFitTable(   modelContainer, dataset, options, ready)
   .efaResidualTable(        modelContainer,          options, ready)
-  .parallelAnalysisTable(   modelContainer, dataset, options, ready, name = "Factor")
-  .efaScreePlot(            modelContainer, dataset, options, ready)
   .efaPathDiagram(          modelContainer, dataset, options, ready)
 
   # data saving
@@ -59,8 +55,8 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
   } else {
     modelContainer <- createJaspContainer()
     modelContainer$dependOn(c("rotationMethod", "orthogonalSelector", "obliqueSelector", "variables", "variables.types",
-                              "factorCountMethod", "eigenvaluesAbove", "manualNumberOfFactors", "naAction",
-                              "baseDecompositionOn", "factoringMethod", "parallelAnalysisMethod", "dataType", "sampleSize"))
+                              "manualNumberOfFactors", "naAction",
+                              "baseDecompositionOn", "factoringMethod", "dataType", "sampleSize"))
     jaspResults[["modelContainer"]] <- modelContainer
   }
 
@@ -142,7 +138,7 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
   efaResult <- try(
     psych::fa(
       r        = dtUse,
-      nfactors = .efaGetNComp(dataset, options, modelContainer),
+      nfactors = options[["manualNumberOfFactors"]],
       rotate   = ifelse(options$rotationMethod == "orthogonal", options$orthogonalSelector, options$obliqueSelector),
       scores   = TRUE,
       covar    = options$baseDecompositionOn == "covarianceMatrix",
@@ -170,61 +166,6 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
 
 }
 
-.efaGetNComp <- function(dataset, options, modelContainer) {
-
-  # manual method
-  if (options$factorCountMethod == "manual") return(options$manualNumberOfFactors)
-
-  # parallel analysis method (do always)
-  if (any(options$variables.types %in% c("ordinal", "nominal")) && options[["baseDecompositionOn"]] == "polyTetrachoricCorrelationMatrix") {
-    polyTetraCor <- .efaComputePolyCorMatrix(modelContainer, dataset, options)
-    if (is.null(polyTetraCor)) return(1)
-
-    .setSeedJASP(options)
-    parallelResult <- try(psych::fa.parallel(polyTetraCor,
-                                             plot = FALSE,
-                                             fa = ifelse(options[["parallelAnalysisMethod"]] == "principalComponentBased",
-                                                         "pc", "fa"),
-                                             n.obs = nrow(dataset)))
-  } else {
-    .setSeedJASP(options)
-
-    parallelResult <- try(psych::fa.parallel(dataset, plot = FALSE,
-                                             fa = ifelse(options[["parallelAnalysisMethod"]] == "principalComponentBased",
-                                                         "pc", "fa")))
-  }
-
-  if (options$factorCountMethod == "parallelAnalysis") {
-    if (isTryError(parallelResult)) {
-      errmsg <- gettextf("Parallel analysis failed. Internal error message: %s", .extractErrorMessage(parallelResult))
-      modelContainer$setError(errmsg)
-      return(1)
-    } else {
-      return(ifelse(options[["parallelAnalysisMethod"]] == "principalComponentBased",
-                    max(1, parallelResult$ncomp), max(1, parallelResult$nfact)))
-
-    }
-  }
-
-  # eigenValue method
-  if (options$factorCountMethod == "eigenvalues") {
-    if (!isTryError(parallelResult)) {
-      ncomp <- sum(parallelResult$pc.values > options$eigenvaluesAbove)
-    } else {
-      ncomp  <- 1
-    }
-    # I can use stop() because it's caught by the try and the message is put on
-    # on the modelcontainer.
-    if (ncomp == 0)
-      .quitAnalysis( gettextf("No factors with an eigenvalue > %1$s. Maximum observed eigenvalue equals %2$s.",
-                              options$eigenvaluesAbove, round(max(parallelResult$fa.values), 3)))
-    return(ncomp)
-  }
-
-
-}
-
-
 # Output functions ----
 .efaKMOtest <- function(modelContainer, dataset, options, ready) {
 
@@ -645,169 +586,6 @@ exploratoryFactorAnalysisInternal <- function(jaspResults, dataset, options, ...
 }
 
 
-.parallelAnalysisTable <- function(modelContainer, dataset, options, ready, name = "Factor") {
-  if (!options[["parallelAnalysisTable"]] || !is.null(modelContainer[["parallelTable"]])) return()
-
-  if (!ready || modelContainer$getError()) return()
-
-  if (any(options$variables.types %in% c("ordinal", "nominal")) && options[["baseDecompositionOn"]] == "polyTetrachoricCorrelationMatrix") {
-    polyTetraCor <- .efaComputePolyCorMatrix(modelContainer, dataset, options)
-    if (is.null(polyTetraCor)) return()
-    .setSeedJASP(options)
-
-    parallelResult <- try(psych::fa.parallel(polyTetraCor,
-                                             plot = FALSE,
-                                             fa = ifelse(options[["parallelAnalysisTableMethod"]] == "principalComponentBased",
-                                                         "pc", "fa"),
-                                             n.obs = nrow(dataset)))
-  } else {
-    .setSeedJASP(options)
-
-    parallelResult <- try(psych::fa.parallel(dataset, plot = FALSE,
-                                             fa = ifelse(options[["parallelAnalysisTableMethod"]] == "principalComponentBased",
-                                                         "pc", "fa")))
-  }
-
-  if (options$parallelAnalysisTableMethod == "principalComponentBased") {
-    eigTitle <- gettext("Real data component eigenvalues")
-    rowsName <- gettext(name)
-    RealDataEigen <- parallelResult$pc.values
-    ResampledEigen <- parallelResult$pc.sim
-    footnote <- gettextf("'*' = %s should be retained. Results from PC-based parallel analysis.", name)
-  } else { # parallelAnalysisMethod is FA
-    eigTitle <- gettext("Real data factor eigenvalues")
-    rowsName <- gettext(name)
-    RealDataEigen <- parallelResult$fa.values
-    ResampledEigen <- parallelResult$fa.sim
-    footnote <- gettextf("'*' = %s should be retained. Results from FA-based parallel analysis.", name)
-  }
-
-  parallelTable <- createJaspTable(gettext("Parallel Analysis"))
-  parallelTable$dependOn(c("parallelAnalysisTable", "parallelAnalysisTableMethod"))
-  parallelTable$addColumnInfo(name = "col", title = "", type = "string")
-
-  parallelTable$addColumnInfo(name = "val1", title = eigTitle, type = "number", format = "dp:3")
-  parallelTable$addColumnInfo(name = "val2", title = gettext("Simulated data mean eigenvalues"), type = "number", format = "dp:3")
-  parallelTable$position <- 6
-  modelContainer[["parallelTable"]] <- parallelTable
-
-
-  efaResult <- modelContainer[["model"]][["object"]]
-  PAs <- zapsmall(as.matrix(data.frame(RealDataEigen, ResampledEigen), fix.empty.names = FALSE))
-  forasterisk <- data.frame(applyforasterisk = RealDataEigen - ResampledEigen)
-  dims <- nrow(PAs)
-
-  firstcol <- data.frame()
-  for (i in 1:dims) {
-    if (forasterisk$applyforasterisk[i] > 0) {
-      firstcol <- rbind(firstcol, paste(rowsName, " ", i,"*", sep = ""))
-    } else {
-      firstcol <- rbind(firstcol, paste(rowsName, i))
-    }
-  }
-
-  parallelTable[["col"]] <- firstcol[[1]]
-  parallelTable[["val1"]] <- c(RealDataEigen)
-  parallelTable[["val2"]] <- c(ResampledEigen)
-
-  parallelTable$addFootnote(message = footnote)
-}
-
-
-.efaScreePlot <- function(modelContainer, dataset, options, ready) {
-  if (!options[["screePlot"]] || !is.null(modelContainer[["scree"]])) return()
-
-  scree <- createJaspPlot(title = "Scree plot", width = 480, height = 320)
-  scree$dependOn(c("screePlot", "screePlotParallelAnalysisResults", "parallelAnalysisMethod"))
-  scree$position <- 8
-  modelContainer[["scree"]] <- scree
-
-  if (!ready || modelContainer$getError()) return()
-
-  n_col <- ncol(dataset)
-
-  if (options[["screePlotParallelAnalysisResults"]]) {
-
-    if (any(options$variables.types %in% c("ordinal", "nominal")) && options[["baseDecompositionOn"]] == "polyTetrachoricCorrelationMatrix") {
-      polyTetraCor <- .efaComputePolyCorMatrix(modelContainer, dataset, options)
-      if (is.null(polyTetraCor)) return()
-      .setSeedJASP(options)
-
-      parallelResult <- try(psych::fa.parallel(polyTetraCor,
-                                   plot = FALSE,
-                                   fa = ifelse(options[["parallelAnalysisTableMethod"]] == "principalComponentBased",
-                                               "pc", "fa"),
-                                   n.obs = nrow(dataset)))
-    } else {
-      .setSeedJASP(options)
-
-      parallelResult <- try(psych::fa.parallel(dataset, plot = FALSE,
-                                               fa = ifelse(options[["parallelAnalysisTableMethod"]] == "principalComponentBased",
-                                                           "pc", "fa")))
-    }
-    if (isTryError(parallelResult)) {
-      errmsg <- gettextf("Screeplot not available. \nInternal error message: %s", .extractErrorMessage(parallelResult))
-      scree$setError(errmsg)
-      # scree$setError(.decodeVarsInMessage(names(dataset), errmsg))
-      return()
-    }
-
-    if (options$parallelAnalysisTableMethod == "factorBased") {
-      evs <- c(parallelResult$fa.values, parallelResult$fa.sim)
-    } else { # in all other cases we use the initial eigenvalues for the plot, aka the pca ones
-      evs <- c(parallelResult$pc.values, parallelResult$pc.sim)
-    }
-    tp <- rep(c(gettext("Data"), gettext("Simulated data from parallel analysis")), each = n_col)
-
-  } else { # do not display parallel analysis
-    evs <- eigen(cor(dataset, use = "pairwise.complete.obs"), only.values = TRUE)$values
-    tp <- rep(gettext("Data"), each = n_col)
-  }
-
-  df <- data.frame(
-    id   = rep(seq_len(n_col), 2),
-    ev   = evs,
-    type = tp
-  )
-
-  # basic scree plot
-  plt <-
-    ggplot2::ggplot(df, ggplot2::aes(x = id, y = ev, linetype = type, shape = type)) +
-    ggplot2::geom_line(na.rm = TRUE) +
-    ggplot2::labs(x = gettext("Factor"), y = gettext("Eigenvalue")) +
-    ggplot2::geom_hline(yintercept = options$eigenvaluesAbove)
-
-
-  # dynamic function for point size:
-  # the plot looks good with size 3 when there are 10 points (3 + log(10) - log(10) = 3)
-  # with more points, the size will become logarithmically smaller until a minimum of
-  # 3 + log(10) - log(200) = 0.004267726
-  # with fewer points, they become bigger to a maximum of 3 + log(10) - log(2) = 4.609438
-  pointsize <- 3 + log(10) - log(n_col)
-  if (pointsize > 0) {
-    plt <- plt + ggplot2::geom_point(na.rm = TRUE, size = max(0, 3 + log(10) - log(n_col)))
-  }
-
-  # add axis lines and better breaks
-  plt <- plt +
-    jaspGraphs::geom_rangeframe() +
-    jaspGraphs::themeJaspRaw() +
-    ggplot2::scale_x_continuous(breaks = seq(1:n_col))
-
-  # theming with special legend thingy
-  plt <- plt +
-    ggplot2::theme(
-      legend.position      = c(0.99, 0.95),
-      legend.justification = c(1, 1),
-      legend.text          = ggplot2::element_text(size = 12.5),
-      legend.title         = ggplot2::element_blank(),
-      legend.key.size      = ggplot2::unit(18, "pt")
-    )
-
-  scree$plotObject <- plt
-  modelContainer[["scree"]] <- scree
-}
-
 .efaPathDiagram <- function(modelContainer, dataset, options, ready){
   if (!options[["pathDiagram"]] || !is.null(modelContainer[["path"]])) return()
 

diff --git a/R/exploratoryfactoranalysisWrapper.R b/R/exploratoryfactoranalysisWrapper.R
@@ -68,8 +68,6 @@
 #'   \item \code{"variables"}
 #'   \item \code{"size"}
 #' }
-#' @param parallelAnalysisTable, If this option is selected, a table will be generated exhibiting a detailed output of the parallel analysis. Can be based on principal component eigenvalues (PC) or factor eigenvalues (FA). The seed is taken from the parallel analysis for determining the number of components/factors above.
-#'    Defaults to \code{FALSE}.
 #' @param pathDiagram, By selecting this option, a visual representation of the direction and strength of the relation between the variable and factor will be displayed.
 #'    Defaults to \code{FALSE}.
 #' @param residualMatrix, Displays a table containing the residual variances and correlations.
@@ -79,10 +77,6 @@
 #'   \item \code{"oblique"}: This method produces components that allow for correlation between the components. This method is selected by default. Several possibilities are available. The default is promax.
 #'   \item \code{"orthogonal"}: This method produces components that are uncorrelated. For this method, there are several possibilities that can be selected.
 #' }
-#' @param screePlot, When selecting this option, a scree plot will be displayed. The scree plot provides information on how much variance in the data, indicated by the eigenvalue, is explained by each factor. A scree plot can be used to decide how many factors should be selected.
-#'    Defaults to \code{FALSE}.
-#' @param screePlotParallelAnalysisResults, Display the results of the parallel analysis in the scree plot. The parallel analysis will be based on PC or FA as defined by the option for the parallel analysis table.
-#'    Defaults to \code{TRUE}.
 #' @param variables, In this box, the variables to perform the analysis on are selected
 exploratoryFactorAnalysis <- function(
           data = NULL,
@@ -93,9 +87,7 @@ exploratoryFactorAnalysis <- function(
           bartlettTest = FALSE,
           baseDecompositionOn = "correlationMatrix",
           dataType = "raw",
-          eigenvaluesAbove = 1,
           factorCorrelations = FALSE,
-          factorCountMethod = "parallelAnalysis",
           factorStructure = FALSE,
           factoringMethod = "minimumResidual",
           fitIndices = FALSE,
@@ -107,19 +99,12 @@ exploratoryFactorAnalysis <- function(
           obliqueSelector = "promax",
           orderLoadingsBy = "size",
           orthogonalSelector = "none",
-          parallelAnalysisMethod = "principalComponentBased",
-          parallelAnalysisTable = FALSE,
-          parallelAnalysisTableMethod = "principalComponentBased",
           pathDiagram = FALSE,
           plotHeight = 320,
           plotWidth = 480,
           residualMatrix = FALSE,
           rotationMethod = "oblique",
           sampleSize = 200,
-          screePlot = FALSE,
-          screePlotParallelAnalysisResults = TRUE,
-          seed = 1,
-          setSeed = FALSE,
           variables = list(types = list(), value = list())) {
 
    defaultArgCalls <- formals(jaspFactor::exploratoryFactorAnalysis)