From 55fcab52f18d8271426bae7d34c3aacfc3fc14ed Mon Sep 17 00:00:00 2001
From: andreassoteriadesmoj
 <91131120+andreassoteriadesmoj@users.noreply.github.com>
Date: Mon, 26 Jan 2026 08:45:54 +0000
Subject: [PATCH 1/4] code: Function `EFA_POOLED()`

---
 R/EFA_POOLED.R | 291 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 291 insertions(+)
 create mode 100644 R/EFA_POOLED.R

diff --git a/R/EFA_POOLED.R b/R/EFA_POOLED.R
new file mode 100644
index 0000000..4d998d6
--- /dev/null
+++ b/R/EFA_POOLED.R
@@ -0,0 +1,291 @@
+#' Exploratory factor analysis on multiple data imputations
+#'
+#' @author Andreas Soteriades
+#'
+#' @description
+#' Consider a dataset with missing values that has been imputed \eqn{m} times.
+#' This function runs \code{\link{EFA}} on each imputed dataset, averages the
+#' estimated parameters (loadings) and calculates confidence intervals (CIs) on
+#' the pooled (mean) loadings.
+#'
+#' @details
+#' Multiple imputation (van Buuren, 2018) was developed by Rubin (1976). The
+#' idea stems from the fact that that imputing one value (single imputation) for
+#' the missing value could not be correct in general. Instead, impute the data
+#' \eqn{m} times, i.e. generate multiple versions of the data table, where each
+#' version has been imputed. Then, fit \eqn{m} models (e.g. \eqn{m} regressions
+#' or \eqn{m} EFA models) and average (pool) the \eqn{m} estimated parameters
+#' (e.g. the \eqn{m} regression coefficients of a given predictor or the \eqn{m}
+#' loadings of a given variable and factor). This average is your estimate.
+#'
+#' In \code{\link{EFA_POOLED}}, the pooled loading \eqn{\bar{\lambda_{ij}}}for
+#' the \eqn{i}th variable and \eqn{j}th factor is calculated as
+#' \eqn{\frac{1}{m}\sum_{k=1}^{m}{\lambda_{ij}^k}}, where \eqn{m} is the number
+#' of imputations
+#'
+#' Similarly, for the \eqn{i}th variable, the communality \code{h2} returned by
+#' \code{\link{EFA_POOLED}} is calculated from the pooled loadings as
+#' \eqn{\sum_{j=1}^{p}{\bar{\lambda_{ij}}^2}}, where \eqn{p} is the number of
+#' factors.
+#'
+#' @source Rubin, D.B. (1976). Inference and Missing Data. Biometrika 63 (3):
+#'   581–90.
+#' @source van Buuren, S. (2018). Flexible Imputation of Missing Data. Second
+#'   Edition. CRC Press. \url{https://stefvanbuuren.name/fimd/}
+#'
+#' @param data_list A list of length \eqn{m}, where \eqn{m} is the number of
+#' imputations. Each list element \eqn{k} is a data frame or matrix of raw data
+#' or a matrix with correlations for data imputation \eqn{k}. See argument
+#' \code{x} in \code{\link{EFA}}.
+#' @param p Numeric. One minus the confidence level of the CIs. Defaults to 0.05
+#' for 95\% CIs.
+#' @param ... Additional arguments passed to \code{\link{EFA}}.
+#'
+#' @return A list of class EFA containing (a subset of) the following:
+#'
+#' \item{h2}{Final communality estimates from the unrotated solution. This is
+#' based on the squared pooled unrotated loadings, \code{unrot_loadings}.}
+#' \item{unrot_loadings}{Loading matrix containing the final, pooled unrotated
+#' loadings.}
+#' \item{vars_accounted}{Matrix of explained variances and sums of squared
+#' loadings. Based on the pooled unrotated loadings.}
+#' \item{fit_indices}{For ML and ULS: Pooled fit indices derived from the unrotated
+#' factor loadings: Chi Square, including significance level, degrees of freedom
+#' (df), Comparative Fit Index (CFI), Root Mean Square Error of Approximation
+#' (RMSEA), including its 90\% confidence interval, Akaike Information Criterion
+#' (AIC), Bayesian Information Criterion (BIC), and the common part accounted
+#' for (CAF) index as proposed by Lorenzo-Seva, Timmerman, & Kiers (2011).
+#' For PAF, only the CAF and dfs are returned.}
+#' \item{rot_loadings}{Loading matrix containing the final pooled rotated
+#' loadings (pattern matrix).}
+#' \item{Phi}{The pooled factor intercorrelations (only for oblique rotations).}
+#' \item{Structure}{The pooled structure matrix (only for oblique rotations).}
+#' \item{vars_accounted_rot}{Matrix of explained variances and sums of squared
+#' loadings. Based on pooled rotated loadings and, for oblique rotations, the
+#' pooled factor intercorrelations.}
+#' \item{settings}{A list of the settings used.}
+#' \item{ci_unrot_loadings}{CIs for the pooled unrotated loadings.}
+#' \item{ci_rot_loadings}{CIs for the pooled rotated loadings.}
+#' \item{ci_structure_loadings}{CIs for the pooled structure loadings (only for
+#' oblique rotations).}
+#' \item{ci_phis}{CIs for the pooled factor intercorrelations (only for oblique
+#' rotations).}
+#' \item{fits}{List of length \eqn{m} with the results from the \code{\link{EFA}}
+#' runs for each imputed dataset.}
+#'
+#' @export
+#'
+#' @examples
+#' dataset <- psych::bfi[1:500, 1:25] # Create dataset
+#' apply(dataset, 2, function(x) sum(is.na(x))) # Any NAs in the columns?
+#'
+#' # Function to impute columns
+#' impute_vector <- function(x, func = mean) {
+#'   replace(
+#'     x,
+#'     is.na(x),
+#'     func(x, na.rm = TRUE)
+#'   )
+#' }
+#'
+#' # Impute data using mean
+#' dataset_imp_mean <- apply(
+#'   dataset,
+#'   2,
+#'   impute_vector,
+#'   func = mean
+#' )
+#'
+#' # Impute data using mean
+#' dataset_imp_median <- apply(
+#'   dataset,
+#'   2,
+#'   impute_vector,
+#'   func = median
+#' )
+#'
+#' # List of imputed datasets, one using the mean, the other using the median
+#' data_list <- list(
+#'   dataset_imp_mean,
+#'   dataset_imp_median
+#' )
+#'
+#' # Run EFA on each imputed dataset, pool the results
+#' efa_pooled <- EFA_POOLED(data_list, p = 0.05, n_factors = 3,
+#'                          type = "EFAtools", method = "ML", rotation = "oblimin")
+#'
+EFA_POOLED <- function(data_list, p = 0.05, ...) {
+  efa_args <- list(...) # Arguments of EFA()
+
+  checkmate::assert_list(data_list)
+  lapply(
+    data_list,
+    checkmate::assert_multi_class,
+    c('matrix', 'data frame')
+  )
+  checkmate::assert_number(p, na.ok = FALSE)
+
+  # Fit EFA for each dataset
+  fits <- lapply(
+    data_list,
+    function(data_list_subset, ...) {
+      # data_list_subset is the i-th element of data_list, and is a data frame
+      # or matrix of raw data or matrix with correlations, i.e. it is argument
+      # x of EFA(). In each lapply() repetition, we modify x to be the i-th
+      # element of data_list. We then run EFA() for this x
+      efa_args$x <- data_list_subset
+      do.call(EFAtools::EFA, efa_args)
+    }
+  )
+
+  # List of structure loadings matrices
+  structure_loadings <- .extract_list_object(fits, 'Structure')
+
+  # For orthogonal rotations, structure_loadings will be empty
+  check_empty <- structure_loadings %>%
+    Filter(Negate(is.null), .) %>%
+    rlang::is_empty()
+
+  # List of unrotated loadings matrices
+  unrot_loadings <- fits %>%
+    .extract_list_object('unrot_loadings') %>%
+    lapply(.change_class, 'matrix')
+
+  # List of rotated loadings matrices
+  rot_loadings <- fits %>%
+    .extract_list_object('rot_loadings') %>%
+    lapply(.change_class, 'matrix')
+
+  # Target-rotate the rotated loadings matrices. This is necessary for averaging
+  # parameters estimated using different imputations
+  target_rotations <- lapply(
+    rot_loadings[-1], # First one is the target- redundant to self-target-rotate
+    psych::target.rot,
+    keys = rot_loadings[[1]]
+  )
+
+  # List of rotated loadings matrices, now target-rotated
+  rot_loadings[-1] <- lapply(
+    target_rotations,
+    function(x) {
+      .change_class(x$loadings, 'matrix')
+    }
+  )
+
+  # Factor correlations
+  if (!check_empty) {
+    phis <- .extract_list_object(fits, 'Phi')
+  }
+
+  # Means of estimated parameters
+  mean_unrot_loadings <- .stat_over_list(unrot_loadings, mean)
+  mean_rot_loadings <- .stat_over_list(rot_loadings, mean)
+  if (!check_empty) {
+    mean_structure_loadings <- .stat_over_list(structure_loadings, mean)
+    mean_phis <- .stat_over_list(phis, mean)
+  }
+
+  # Standard deviations of estimated parameters
+  sd_unrot_loadings <- .stat_over_list(unrot_loadings, sd)
+  sd_rot_loadings <- .stat_over_list(rot_loadings, sd)
+  if (!check_empty) {
+    sd_structure_loadings <- .stat_over_list(structure_loadings, sd)
+    sd_phis <- .stat_over_list(phis, sd)
+  }
+
+  # CIs of estimated parameters
+  n_perm <- length(fits)
+  ci_unrot_loadings <- .calc_cis(
+    mean_unrot_loadings,
+    sd_unrot_loadings,
+    p,
+    n = n_perm
+  )
+
+  ci_rot_loadings <- .calc_cis(
+    mean_rot_loadings,
+    sd_rot_loadings,
+    p,
+    n = n_perm
+  )
+
+  if (!check_empty) {
+    ci_structure_loadings <- .calc_cis(
+      mean_structure_loadings,
+      sd_structure_loadings,
+      p,
+      n = n_perm
+    )
+
+    ci_phis <- .calc_cis(
+      mean_phis,
+      sd_phis,
+      p,
+      n = n_perm
+    )
+  }
+
+  # Communalities from mean unrotated loadings
+  mean_h2 <- rowSums(mean_unrot_loadings ^ 2)
+
+  # Variances, mean unrotated loadings
+  mean_vars_accounted <- EFAtools:::.compute_vars(
+    L_unrot = mean_unrot_loadings,
+    L_rot = mean_unrot_loadings,
+    Phi = NULL
+  )
+
+  # Variances, mean rotated loadings
+  mean_vars_accounted_rot <- EFAtools:::.compute_vars(
+    L_unrot = mean_rot_loadings,
+    L_rot = mean_rot_loadings,
+    Phi = mean_phis
+  )
+
+  mean_fit_indices <- fits %>%
+    .extract_list_object('fit_indices') %>%
+    lapply(unlist) %>%
+    simplify2array() %>%
+    rowMeans() %>%
+    as.list()
+
+  # Ensure loadings tables class is consistent with EFA() output
+  mean_unrot_loadings <- .change_class(mean_unrot_loadings, "LOADINGS")
+  mean_rot_loadings <- .change_class(mean_rot_loadings, "LOADINGS")
+
+  if (!check_empty) {
+    results <- list(
+      h2 = mean_h2,
+      unrot_loadings = mean_unrot_loadings,
+      vars_accounted = mean_vars_accounted,
+      fit_indices = mean_fit_indices,
+      rot_loadings = mean_rot_loadings,
+      Phi = mean_phis,
+      Structure = mean_structure_loadings,
+      vars_accounted_rot = mean_vars_accounted_rot,
+      settings = fits[[1]]$settings,
+      ci_unrot_loadings = ci_unrot_loadings,
+      ci_rot_loadings = ci_rot_loadings,
+      ci_structure_loadings = ci_structure_loadings,
+      ci_phis = ci_phis,
+      fits = fits
+    )
+  } else {
+    results <- list(
+      h2 = mean_h2,
+      unrot_loadings = mean_unrot_loadings,
+      vars_accounted = mean_vars_accounted,
+      fit_indices = mean_fit_indices,
+      rot_loadings = mean_rot_loadings,
+      vars_accounted_rot = mean_vars_accounted_rot,
+      settings = fits[[1]]$settings,
+      ci_unrot_loadings = ci_unrot_loadings,
+      ci_rot_loadings = ci_rot_loadings,
+      fits = fits
+    )
+  }
+
+  # class(results) <- 'EFA'
+
+  return(results)
+}

From 2e022b8c115ef4abcf96f3efc8a518097d4c79a2 Mon Sep 17 00:00:00 2001
From: andreassoteriadesmoj
 <91131120+andreassoteriadesmoj@users.noreply.github.com>
Date: Mon, 26 Jan 2026 08:47:55 +0000
Subject: [PATCH 2/4] code: Helper functions for internal use by `EFA_POOLED()`

---
 R/helper.R | 111 ++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 110 insertions(+), 1 deletion(-)

diff --git a/R/helper.R b/R/helper.R
index cdfbb09..30d2ec3 100644
--- a/R/helper.R
+++ b/R/helper.R
@@ -1,3 +1,113 @@
+#' Extract a list object by its name
+#'
+#' Consider a list of named sub-lists. This function extracts, for each sub-list,
+#' the sub-list element that is specified by the user. This function is useful
+#' for extracting results from \code{\link{EFA}} for each permutation run in
+#' \code{\link{EFA_POOLED}}.
+#'
+#' @param alist A list of sub-lists, typically a list of \eqn{m} objects of class
+#' \code{"EFA"}, where \eqn{m} is the number of imputations passed to
+#' \code{\link{EFA_POOLED}}.
+#' @param object String of length 1. The name of the object to extract e.g.
+#' \code{"h2"} or \code{"vars_accounted"}.
+#'
+#' @return A list of length \eqn{m}, with each element containing the extracted
+#' \code{object} for the \eqn{k}th element (\eqn{k = 1,..., m}).
+.extract_list_object <- function(alist, object) {
+  lapply(
+    alist,
+    function(x) {
+      x[[object]]
+    }
+  )
+}
+
+#' Calculate statistics for a list of matrices
+#'
+#' Given a list of matrices, this function calculates user-supplied statistics
+#' (e.g. mean, median) over the matrices. This function is useful for averaging
+#' results from \code{\link{EFA}} (e.g. loadings) over all permutation runs in
+#' \code{\link{EFA_POOLED}}.
+#'
+#' @param alist A list of sub-lists, typically a list of \eqn{m} matrices from
+#' \code{"EFA"}, where \eqn{m} is the number of imputations passed to
+#' \code{\link{EFA_POOLED}}.
+#' @param stat A function, e.g. \code{mean} or \code{sd}.
+#'
+#' @return A matrix with the aggregated results
+#' @examples
+#' alist <- list(
+#'   matrix(1:10, 2),
+#'   matrix(11:20, 2)
+#'   )
+#'
+#' # Have a look at list
+#' alist
+#'
+#' # Mean matrix
+#' .stat_over_list(alist, mean)
+#'
+#' # Same as averaging the two matrices in the list
+#' identical(
+#'   .stat_over_list(alist, mean),
+#'   (alist[[1]] + alist[[2]]) / 2
+#' )
+.stat_over_list <- function(alist, stat) {
+  apply(
+    simplify2array(alist),
+    seq_along(alist),
+    stat
+  )
+}
+
+#' Covert a \code{"LOADINGS"} table to matrix or a matrix to \code{"LOADINGS"}
+#'
+#' The loadings tables returned by \code{\link{EFA}} are of class
+#' \code{"LOADINGS"}, which prevents applying functions on them. This function
+#' allows to change their class to \code{"matrix"}, and to change back to
+#' \code{"LOADINGS"} when done.
+#'
+#' @param x A table of class \code{"matrix"} or \code{"LOADINGS"}.
+#' @param cl A string with the class to change the table to. Should be
+#' \code{"LOADINGS"} or \code{"matrix"}.
+#'
+#' @return A table with the loadings, of class either \code{"LOADINGS"} or
+#' \code{"matrix"}.
+.change_class <- function(x, cl = 'matrix') {
+  class(x) <- cl
+  return(x)
+}
+
+#' Confidence intervals around mean
+#'
+#' This function is used internally by \code{\link{EFA_POOLED}} to calculates
+#' confidence intervals (CIs) around the pooled loadings and pooled interfactor
+#' correlations.
+#'
+#' @param means A matrix with the pooled loadings or pooled interfactor
+#' correlations
+#' @param sds A matrix with the standard deviations for the loadings or
+#' interfactor correlations
+#' @param p Numeric. One minus the confidence level of the CIs. Defaults to 0.05
+#' for 95\% CIs.
+#' @param n Numeric. Typically, the number of permutations \eqn{m}. See
+#' \code{\link{EFA_POOLED}}
+#'
+#' @return A list with the lower and upper CIs.
+.calc_cis <- function(means, sds, p = 0.05, n) {
+  error <- qnorm(1 - p / 2) * sds / sqrt(n)
+
+  lower_ci <- means - error
+  upper_ci <-  means + error
+
+  cis <- list(
+    lower = lower_ci,
+    upper = upper_ci
+  )
+
+  return(cis)
+}
+
 #' Format numbers for print method
 #'
 #' Helper function used in the print method for class LOADINGS and SLLOADINGS.
@@ -1060,4 +1170,3 @@ if(n == 1){
   return(do.call(rbind, t_grid_list))
 }
 
-

From a02df0710483144529275a8a6c1c7a7e4f85fd23 Mon Sep 17 00:00:00 2001
From: andreassoteriadesmoj
 <91131120+andreassoteriadesmoj@users.noreply.github.com>
Date: Mon, 26 Jan 2026 08:48:51 +0000
Subject: [PATCH 3/4] doc: Add new function documentation

`EFA_POOLED()` & new helper functions
---
 man/EFA_POOLED.Rd              | 128 +++++++++++++++++++++++++++++++++
 man/dot-calc_cis.Rd            |  29 ++++++++
 man/dot-change_class.Rd        |  24 +++++++
 man/dot-extract_list_object.Rd |  26 +++++++
 man/dot-stat_over_list.Rd      |  42 +++++++++++
 5 files changed, 249 insertions(+)
 create mode 100644 man/EFA_POOLED.Rd
 create mode 100644 man/dot-calc_cis.Rd
 create mode 100644 man/dot-change_class.Rd
 create mode 100644 man/dot-extract_list_object.Rd
 create mode 100644 man/dot-stat_over_list.Rd

diff --git a/man/EFA_POOLED.Rd b/man/EFA_POOLED.Rd
new file mode 100644
index 0000000..c3c0e36
--- /dev/null
+++ b/man/EFA_POOLED.Rd
@@ -0,0 +1,128 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/EFA_POOLED.R
+\name{EFA_POOLED}
+\alias{EFA_POOLED}
+\title{Exploratory factor analysis on multiple data imputations}
+\source{
+Rubin, D.B. (1976). Inference and Missing Data. Biometrika 63 (3):
+  581–90.
+
+van Buuren, S. (2018). Flexible Imputation of Missing Data. Second
+  Edition. CRC Press. \url{https://stefvanbuuren.name/fimd/}
+}
+\usage{
+EFA_POOLED(data_list, p = 0.05, ...)
+}
+\arguments{
+\item{data_list}{A list of length \eqn{m}, where \eqn{m} is the number of
+imputations. Each list element \eqn{k} is a data frame or matrix of raw data
+or a matrix with correlations for data imputation \eqn{k}. See argument
+\code{x} in \code{\link{EFA}}.}
+
+\item{p}{Numeric. One minus the confidence level of the CIs. Defaults to 0.05
+for 95\% CIs.}
+
+\item{...}{Additional arguments passed to \code{\link{EFA}}.}
+}
+\value{
+A list of class EFA containing (a subset of) the following:
+
+\item{h2}{Final communality estimates from the unrotated solution. This is
+based on the squared pooled unrotated loadings, \code{unrot_loadings}.}
+\item{unrot_loadings}{Loading matrix containing the final, pooled unrotated
+loadings.}
+\item{vars_accounted}{Matrix of explained variances and sums of squared
+loadings. Based on the pooled unrotated loadings.}
+\item{fit_indices}{For ML and ULS: Pooled fit indices derived from the unrotated
+factor loadings: Chi Square, including significance level, degrees of freedom
+(df), Comparative Fit Index (CFI), Root Mean Square Error of Approximation
+(RMSEA), including its 90\% confidence interval, Akaike Information Criterion
+(AIC), Bayesian Information Criterion (BIC), and the common part accounted
+for (CAF) index as proposed by Lorenzo-Seva, Timmerman, & Kiers (2011).
+For PAF, only the CAF and dfs are returned.}
+\item{rot_loadings}{Loading matrix containing the final pooled rotated
+loadings (pattern matrix).}
+\item{Phi}{The pooled factor intercorrelations (only for oblique rotations).}
+\item{Structure}{The pooled structure matrix (only for oblique rotations).}
+\item{vars_accounted_rot}{Matrix of explained variances and sums of squared
+loadings. Based on pooled rotated loadings and, for oblique rotations, the
+pooled factor intercorrelations.}
+\item{settings}{A list of the settings used.}
+\item{ci_unrot_loadings}{CIs for the pooled unrotated loadings.}
+\item{ci_rot_loadings}{CIs for the pooled rotated loadings.}
+\item{ci_structure_loadings}{CIs for the pooled structure loadings (only for
+oblique rotations).}
+\item{ci_phis}{CIs for the pooled factor intercorrelations (only for oblique
+rotations).}
+\item{fits}{List of length \eqn{m} with the results from the \code{\link{EFA}}
+runs for each imputed dataset.}
+}
+\description{
+Consider a dataset with missing values that has been imputed \eqn{m} times.
+This function runs \code{\link{EFA}} on each imputed dataset, averages the
+estimated parameters (loadings) and calculates confidence intervals (CIs) on
+the pooled (mean) loadings.
+}
+\details{
+Multiple imputation (van Buuren, 2018) was developed by Rubin (1976). The
+idea stems from the fact that that imputing one value (single imputation) for
+the missing value could not be correct in general. Instead, impute the data
+\eqn{m} times, i.e. generate multiple versions of the data table, where each
+version has been imputed. Then, fit \eqn{m} models (e.g. \eqn{m} regressions
+or \eqn{m} EFA models) and average (pool) the \eqn{m} estimated parameters
+(e.g. the \eqn{m} regression coefficients of a given predictor or the \eqn{m}
+loadings of a given variable and factor). This average is your estimate.
+
+In \code{\link{EFA_POOLED}}, the pooled loading \eqn{\bar{\lambda_{ij}}}for
+the \eqn{i}th variable and \eqn{j}th factor is calculated as
+\eqn{\frac{1}{m}\sum_{k=1}^{m}{\lambda_{ij}^k}}, where \eqn{m} is the number
+of imputations
+
+Similarly, for the \eqn{i}th variable, the communality \code{h2} returned by
+\code{\link{EFA_POOLED}} is calculated from the pooled loadings as
+\eqn{\sum_{j=1}^{p}{\bar{\lambda_{ij}}^2}}, where \eqn{p} is the number of
+factors.
+}
+\examples{
+dataset <- psych::bfi[1:500, 1:25] # Create dataset
+apply(dataset, 2, function(x) sum(is.na(x))) # Any NAs in the columns?
+
+# Function to impute columns
+impute_vector <- function(x, func = mean) {
+  replace(
+    x,
+    is.na(x),
+    func(x, na.rm = TRUE)
+  )
+}
+
+# Impute data using mean
+dataset_imp_mean <- apply(
+  dataset,
+  2,
+  impute_vector,
+  func = mean
+)
+
+# Impute data using mean
+dataset_imp_median <- apply(
+  dataset,
+  2,
+  impute_vector,
+  func = median
+)
+
+# List of imputed datasets, one using the mean, the other using the median
+data_list <- list(
+  dataset_imp_mean,
+  dataset_imp_median
+)
+
+# Run EFA on each imputed dataset, pool the results
+efa_pooled <- EFA_POOLED(data_list, p = 0.05, n_factors = 3,
+                         type = "EFAtools", method = "ML", rotation = "oblimin")
+
+}
+\author{
+Andreas Soteriades
+}
diff --git a/man/dot-calc_cis.Rd b/man/dot-calc_cis.Rd
new file mode 100644
index 0000000..cad6c14
--- /dev/null
+++ b/man/dot-calc_cis.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper.R
+\name{.calc_cis}
+\alias{.calc_cis}
+\title{Confidence intervals around mean}
+\usage{
+.calc_cis(means, sds, p = 0.05, n)
+}
+\arguments{
+\item{means}{A matrix with the pooled loadings or pooled interfactor
+correlations}
+
+\item{sds}{A matrix with the standard deviations for the loadings or
+interfactor correlations}
+
+\item{p}{Numeric. One minus the confidence level of the CIs. Defaults to 0.05
+for 95\% CIs.}
+
+\item{n}{Numeric. Typically, the number of permutations \eqn{m}. See
+\code{\link{EFA_POOLED}}}
+}
+\value{
+A list with the lower and upper CIs.
+}
+\description{
+This function is used internally by \code{\link{EFA_POOLED}} to calculates
+confidence intervals (CIs) around the pooled loadings and pooled interfactor
+correlations.
+}
diff --git a/man/dot-change_class.Rd b/man/dot-change_class.Rd
new file mode 100644
index 0000000..c46c16d
--- /dev/null
+++ b/man/dot-change_class.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper.R
+\name{.change_class}
+\alias{.change_class}
+\title{Covert a \code{"LOADINGS"} table to matrix or a matrix to \code{"LOADINGS"}}
+\usage{
+.change_class(x, cl = "matrix")
+}
+\arguments{
+\item{x}{A table of class \code{"matrix"} or \code{"LOADINGS"}.}
+
+\item{cl}{A string with the class to change the table to. Should be
+\code{"LOADINGS"} or \code{"matrix"}.}
+}
+\value{
+A table with the loadings, of class either \code{"LOADINGS"} or
+\code{"matrix"}.
+}
+\description{
+The loadings tables returned by \code{\link{EFA}} are of class
+\code{"LOADINGS"}, which prevents applying functions on them. This function
+allows to change their class to \code{"matrix"}, and to change back to
+\code{"LOADINGS"} when done.
+}
diff --git a/man/dot-extract_list_object.Rd b/man/dot-extract_list_object.Rd
new file mode 100644
index 0000000..e231471
--- /dev/null
+++ b/man/dot-extract_list_object.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper.R
+\name{.extract_list_object}
+\alias{.extract_list_object}
+\title{Extract a list object by its name}
+\usage{
+.extract_list_object(alist, object)
+}
+\arguments{
+\item{alist}{A list of sub-lists, typically a list of \eqn{m} objects of class
+\code{"EFA"}, where \eqn{m} is the number of imputations passed to
+\code{\link{EFA_POOLED}}.}
+
+\item{object}{String of length 1. The name of the object to extract e.g.
+\code{"h2"} or \code{"vars_accounted"}.}
+}
+\value{
+A list of length \eqn{m}, with each element containing the extracted
+\code{object} for the \eqn{k}th element (\eqn{k = 1,..., m}).
+}
+\description{
+Consider a list of named sub-lists. This function extracts, for each sub-list,
+the sub-list element that is specified by the user. This function is useful
+for extracting results from \code{\link{EFA}} for each permutation run in
+\code{\link{EFA_POOLED}}.
+}
diff --git a/man/dot-stat_over_list.Rd b/man/dot-stat_over_list.Rd
new file mode 100644
index 0000000..f47b7ed
--- /dev/null
+++ b/man/dot-stat_over_list.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper.R
+\name{.stat_over_list}
+\alias{.stat_over_list}
+\title{Calculate statistics for a list of matrices}
+\usage{
+.stat_over_list(alist, stat)
+}
+\arguments{
+\item{alist}{A list of sub-lists, typically a list of \eqn{m} matrices from
+\code{"EFA"}, where \eqn{m} is the number of imputations passed to
+\code{\link{EFA_POOLED}}.}
+
+\item{stat}{A function, e.g. \code{mean} or \code{sd}.}
+}
+\value{
+A matrix with the aggregated results
+}
+\description{
+Given a list of matrices, this function calculates user-supplied statistics
+(e.g. mean, median) over the matrices. This function is useful for averaging
+results from \code{\link{EFA}} (e.g. loadings) over all permutation runs in
+\code{\link{EFA_POOLED}}.
+}
+\examples{
+alist <- list(
+  matrix(1:10, 2),
+  matrix(11:20, 2)
+  )
+
+# Have a look at list
+alist
+
+# Mean matrix
+.stat_over_list(alist, mean)
+
+# Same as averaging the two matrices in the list
+identical(
+  .stat_over_list(alist, mean),
+  (alist[[1]] + alist[[2]]) / 2
+)
+}

From cb076dbdc733d17510f7303a7b8113bcdfcb316c Mon Sep 17 00:00:00 2001
From: andreassoteriadesmoj
 <91131120+andreassoteriadesmoj@users.noreply.github.com>
Date: Mon, 26 Jan 2026 08:49:04 +0000
Subject: [PATCH 4/4] Update NAMESPACE

---
 NAMESPACE | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NAMESPACE b/NAMESPACE
index d2c6ce1..249eb9f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -32,6 +32,7 @@ export(CD)
 export(COMPARE)
 export(EFA)
 export(EFA_AVERAGE)
+export(EFA_POOLED)
 export(EKC)
 export(FACTOR_SCORES)
 export(HULL)