Nothing
#' Number of components/factors to retain in PCA/FA
#'
#' This function runs many existing procedures for determining how many factors
#' to retain/extract from factor analysis (FA) or dimension reduction (PCA). It
#' returns the number of factors based on the maximum consensus between methods.
#' In case of ties, it will keep the simplest model and select the solution
#' with the fewer factors.
#'
#' @param x A data frame.
#' @param type Can be `"FA"` or `"PCA"`, depending on what you want to do.
#' @param rotation Only used for VSS (Very Simple Structure criterion, see
#' [psych::VSS()]). The rotation to apply. Can be `"none"`, `"varimax"`,
#' `"quartimax"`, `"bentlerT"`, `"equamax"`, `"varimin"`, `"geominT"` and
#' `"bifactor"` for orthogonal rotations, and `"promax"`, `"oblimin"`,
#' `"simplimax"`, `"bentlerQ"`, `"geominQ"`, `"biquartimin"` and `"cluster"`
#' for oblique transformations.
#' @param algorithm Factoring method used by VSS. Can be `"pa"` for Principal
#' Axis Factor Analysis, `"minres"` for minimum residual (OLS) factoring,
#' `"mle"` for Maximum Likelihood FA and `"pc"` for Principal Components.
#' `"default"` will select `"minres"` if `type = "FA"` and `"pc"` if
#' `type = "PCA"`.
#' @param package Package from which respective methods are used. Can be
#' `"all"` or a vector containing `"nFactors"`, `"psych"`, `"PCDimension"`,
#' `"fit"` or `"EGAnet"`. Note that `"fit"` (which actually also relies on the
#' `psych` package) and `"EGAnet"` can be very slow for bigger datasets. Thus,
#' the default is `c("nFactors", "psych")`. You must have the respective
#' packages installed for the methods to be used.
#' @param safe If `TRUE`, the function will run all the procedures in try
#' blocks, and will only return those that work and silently skip the ones
#' that may fail.
#' @param cor An optional correlation matrix that can be used (note that the
#' data must still be passed as the first argument). If `NULL`, will
#' compute it by running `cor()` on the passed data.
#' @param n_max If set to a value (e.g., `10`), will drop from the results all
#' methods that suggest a higher number of components. The interpretation becomes
#' 'from all the methods that suggested a number lower than n_max, the results
#' are ...'.
#' @param ... Arguments passed to or from other methods.
#'
#' @details `n_components()` is actually an alias for `n_factors()`, with
#' different defaults for the function arguments.
#'
#' @note There is also a
#' [`plot()`-method](https://easystats.github.io/see/articles/parameters.html)
#' implemented in the [**see**-package](https://easystats.github.io/see/).
#' `n_components()` is a convenient short-cut for `n_factors(type = "PCA")`.
#'
#' @examplesIf require("PCDimension", quietly = TRUE) && require("nFactors", quietly = TRUE) && require("EGAnet", quietly = TRUE) && require("psych", quietly = TRUE)
#' library(parameters)
#' n_factors(mtcars, type = "PCA")
#'
#' result <- n_factors(mtcars[1:5], type = "FA")
#' as.data.frame(result)
#' summary(result)
#' \donttest{
#' # Setting package = 'all' will increase the number of methods (but is slow)
#' n_factors(mtcars, type = "PCA", package = "all")
#' n_factors(mtcars, type = "FA", algorithm = "mle", package = "all")
#' }
#'
#' @return A data frame.
#'
#' @references
#'
#' - Bartlett, M. S. (1950). Tests of significance in factor analysis.
#' British Journal of statistical psychology, 3(2), 77-85.
#'
#' - Bentler, P. M., & Yuan, K. H. (1996). Test of linear trend in
#' eigenvalues of a covariance matrix with application to data analysis.
#' British Journal of Mathematical and Statistical Psychology, 49(2), 299-312.
#'
#' - Cattell, R. B. (1966). The scree test for the number of factors.
#' Multivariate behavioral research, 1(2), 245-276.
#'
#' - Finch, W. H. (2019). Using Fit Statistic Differences to Determine the
#' Optimal Number of Factors to Retain in an Exploratory Factor Analysis.
#' Educational and Psychological Measurement.
#'
#' - Zoski, K. W., & Jurs, S. (1996). An objective counterpart to the
#' visual scree test for factor analysis: The standard error scree.
#' Educational and Psychological Measurement, 56(3), 443-451.
#'
#' - Zoski, K., & Jurs, S. (1993). Using multiple regression to determine
#' the number of factors to retain in factor analysis. Multiple Linear
#' Regression Viewpoints, 20(1), 5-9.
#'
#' - Nasser, F., Benson, J., & Wisenbaker, J. (2002). The performance of
#' regression-based variations of the visual scree for determining the number
#' of common factors. Educational and psychological measurement, 62(3),
#' 397-419.
#'
#' - Golino, H., Shi, D., Garrido, L. E., Christensen, A. P., Nieto, M.
#' D., Sadana, R., & Thiyagarajan, J. A. (2018). Investigating the performance
#' of Exploratory Graph Analysis and traditional techniques to identify the
#' number of latent factors: A simulation and tutorial.
#'
#' - Golino, H. F., & Epskamp, S. (2017). Exploratory graph analysis: A
#' new approach for estimating the number of dimensions in psychological
#' research. PloS one, 12(6), e0174035.
#'
#' - Revelle, W., & Rocklin, T. (1979). Very simple structure: An
#' alternative procedure for estimating the optimal number of interpretable
#' factors. Multivariate Behavioral Research, 14(4), 403-414.
#'
#' - Velicer, W. F. (1976). Determining the number of components from the
#' matrix of partial correlations. Psychometrika, 41(3), 321-327.
#'
#' @export
n_factors <- function(x,
type = "FA",
rotation = "varimax",
algorithm = "default",
package = c("nFactors", "psych"),
cor = NULL,
safe = TRUE,
n_max = NULL,
...) {
if (all(package == "all")) {
package <- c("nFactors", "EGAnet", "psych", "fit", "pcdimension")
}
# Get number of observations
if (is.data.frame(x)) {
n_obs <- nrow(x)
} else if (is.numeric(x) && !is.null(cor)) {
n_obs <- x
package <- package[!package %in% c("pcdimension", "PCDimension")]
} else if (is.matrix(x) || inherits(x, "easycormatrix")) {
insight::format_error(
"Please input the correlation matrix via the `cor` argument and the number of rows / observations via the first argument." # nolint
)
}
# Get only numeric
numerics <- vapply(x, is.numeric, TRUE)
if (!all(numerics)) {
insight::format_warning(paste0(
"Some variables are not numeric (",
toString(names(x)[!numerics]),
"). Dropping them."
))
}
x <- x[numerics]
# Correlation matrix
if (is.null(cor)) {
cor <- stats::cor(x, use = "pairwise.complete.obs", ...)
}
eigen_values <- eigen(cor)$values
# Smooth matrix if negative eigen values
if (any(eigen_values < 0)) {
insight::check_if_installed("psych")
cor <- psych::cor.smooth(cor, ...)
eigen_values <- eigen(cor)$values
}
# Initialize dataframe
out <- data.frame()
# nFactors -------------------------------------------
if ("nFactors" %in% package) {
insight::check_if_installed("nFactors")
# Model
if (tolower(type) %in% c("fa", "factor", "efa")) {
model <- "factors"
} else {
model <- "components"
}
# Compute all
if (safe) {
out <- rbind(
out,
tryCatch(.n_factors_bartlett(eigen_values, model, n_obs),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
out <- rbind(
out,
tryCatch(.n_factors_bentler(eigen_values, model, n_obs),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
out <- rbind(
out,
tryCatch(.n_factors_cng(eigen_values, model),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
out <- rbind(
out,
tryCatch(.n_factors_mreg(eigen_values, model),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
out <- rbind(
out,
tryCatch(.n_factors_scree(eigen_values, model),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
out <- rbind(
out,
tryCatch(.n_factors_sescree(eigen_values, model),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
} else {
out <- rbind(
out,
.n_factors_bartlett(eigen_values, model, n_obs)
)
out <- rbind(
out,
.n_factors_bentler(eigen_values, model, n_obs)
)
out <- rbind(
out,
.n_factors_cng(eigen_values, model)
)
out <- rbind(
out,
.n_factors_mreg(eigen_values, model)
)
out <- rbind(
out,
.n_factors_scree(eigen_values, model)
)
out <- rbind(
out,
.n_factors_sescree(eigen_values, model)
)
}
}
# EGAnet -------------------------------------------
if ("EGAnet" %in% package) {
insight::check_if_installed("EGAnet")
if (safe) {
out <- rbind(
out,
tryCatch(.n_factors_ega(x, cor, n_obs, eigen_values, type),
# warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
} else {
out <- rbind(
out,
.n_factors_ega(x, cor, n_obs, eigen_values, type)
)
}
}
# psych -------------------------------------------
if ("psych" %in% package) {
insight::check_if_installed("psych")
if (safe) {
out <- rbind(
out,
tryCatch(.n_factors_vss(x, cor, n_obs, type, rotation, algorithm),
# warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
} else {
out <- rbind(
out,
.n_factors_vss(x, cor, n_obs, type, rotation, algorithm)
)
}
}
# fit -------------------------------------------
if ("fit" %in% package) {
insight::check_if_installed("psych")
if (safe) {
out <- rbind(
out,
tryCatch(.n_factors_fit(x, cor, n_obs, type, rotation, algorithm),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
} else {
out <- rbind(
out,
.n_factors_fit(x, cor, n_obs, type, rotation, algorithm)
)
}
}
# pcdimension -------------------------------------------
if ("pcdimension" %in% tolower(package)) {
insight::check_if_installed("PCDimension")
if (safe) {
out <- rbind(
out,
tryCatch(.n_factors_PCDimension(x, type),
warning = function(w) data.frame(),
error = function(e) data.frame()
)
)
} else {
out <- rbind(
out,
.n_factors_PCDimension(x, type)
)
}
}
# OUTPUT ----------------------------------------------
# TODO created weighted composite score
out <- out[!is.na(out$n_Factors), ] # Remove empty methods
out <- out[order(out$n_Factors), ] # Arrange by n factors
row.names(out) <- NULL # Reset row index
if (!is.null(n_max)) {
out <- out[out$n_Factors <= n_max, ]
}
# Add summary
by_factors <- .data_frame(
n_Factors = as.numeric(unique(out$n_Factors)),
n_Methods = as.numeric(by(out, as.factor(out$n_Factors), function(out) n <- nrow(out)))
)
# Add cumulative percentage of variance explained
fa <- factor_analysis(x, cor = cor, n = max(by_factors$n_Factors)) # Get it from our fa:: wrapper (TODO: that's probably not the most efficient)
varex <- attributes(fa)$summary
# Extract number of factors from EFA output (usually MR1, ML1, etc.)
varex$n_Factors <- as.numeric(gsub("[^\\d]+", "", varex$Component, perl = TRUE))
# Merge (and like that filter out empty methods)
by_factors <- merge(by_factors, varex[, c("n_Factors", "Variance_Cumulative")], by = "n_Factors")
attr(out, "Variance_Explained") <- varex # We add all the variance explained (for plotting)
attr(out, "summary") <- by_factors
attr(out, "n") <- min(as.numeric(as.character(
by_factors[by_factors$n_Methods == max(by_factors$n_Methods), "n_Factors"]
)))
class(out) <- c("n_factors", "see_n_factors", class(out))
out
}
#' @rdname n_factors
#' @export
n_components <- function(x,
type = "PCA",
rotation = "varimax",
algorithm = "default",
package = c("nFactors", "psych"),
cor = NULL,
safe = TRUE,
...) {
n_factors(
x,
type = type,
rotation = rotation,
algorithm = algorithm,
package = package,
cor = cor,
safe = safe,
...
)
}
#' @export
print.n_factors <- function(x, ...) {
results <- attributes(x)$summary
# Extract info
max_methods <- max(results$n_Methods)
best_n <- attributes(x)$n
# Extract methods
if ("n_Factors" %in% names(x)) {
type <- "factor"
methods_text <- toString(as.character(x[x$n_Factors == best_n, "Method"]))
} else {
type <- "cluster"
methods_text <- toString(as.character(x[x$n_Clusters == best_n, "Method"]))
}
# Text
msg_text <- paste0(
"The choice of ",
as.character(best_n),
ifelse(type == "factor", " dimensions ", " clusters "),
"is supported by ",
max_methods,
" (",
sprintf("%.2f", max_methods / nrow(x) * 100),
"%) methods out of ",
nrow(x),
" (",
methods_text,
").\n"
)
insight::print_color("# Method Agreement Procedure:\n\n", "blue")
cat(msg_text)
invisible(x)
}
#' @export
summary.n_factors <- function(object, ...) {
attributes(object)$summary
}
#' @export
as.numeric.n_factors <- function(x, ...) {
attributes(x)$n
}
#' @export
as.double.n_factors <- as.numeric.n_factors
#' @export
summary.n_clusters <- summary.n_factors
#' @export
as.numeric.n_clusters <- as.numeric.n_factors
#' @export
as.double.n_clusters <- as.double.n_factors
#' @export
print.n_clusters <- print.n_factors
# Methods -----------------------------------------------------------------
#' Bartlett, Anderson and Lawley Procedures
#' @keywords internal
.n_factors_bartlett <- function(eigen_values = NULL, model = "factors", nobs = NULL) {
nfac <- nFactors::nBartlett(
eigen_values,
N = nobs,
cor = TRUE,
alpha = 0.05,
details = FALSE
)$nFactors
.data_frame(
n_Factors = as.numeric(nfac),
Method = insight::format_capitalize(names(nfac)),
Family = "Barlett"
)
}
#' Bentler and Yuan's Procedure
#' @keywords internal
.n_factors_bentler <- function(eigen_values = NULL, model = "factors", nobs = NULL) {
nfac <- .nBentler(
x = eigen_values,
N = nobs,
model = model,
alpha = 0.05,
details = FALSE
)$nFactors
.data_frame(
n_Factors = as.numeric(nfac),
Method = "Bentler",
Family = "Bentler"
)
}
#' Cattell-Nelson-Gorsuch CNG Indices
#' @keywords internal
.n_factors_cng <- function(eigen_values = NULL, model = "factors") {
if (length(eigen_values) < 6) {
nfac <- NA
} else {
nfac <- nFactors::nCng(x = eigen_values, cor = TRUE, model = model)$nFactors
}
.data_frame(
n_Factors = as.numeric(nfac),
Method = "CNG",
Family = "CNG"
)
}
#' Multiple Regression Procedure
#' @keywords internal
.n_factors_mreg <- function(eigen_values = NULL, model = "factors") {
if (length(eigen_values) < 6) {
nfac <- NA
} else {
nfac <- nFactors::nMreg(x = eigen_values, cor = TRUE, model = model)$nFactors
}
.data_frame(
n_Factors = as.numeric(nfac),
Method = c("beta", "t", "p"),
Family = "Multiple_regression"
)
}
#' Non Graphical Cattell's Scree Test
#' @keywords internal
.n_factors_scree <- function(eigen_values = NULL, model = "factors") {
nfac <- unlist(nFactors::nScree(x = eigen_values, cor = TRUE, model = model)$Components)
.data_frame(
n_Factors = as.numeric(nfac),
Method = c("Optimal coordinates", "Acceleration factor", "Parallel analysis", "Kaiser criterion"),
Family = "Scree"
)
}
#' Standard Error Scree and Coefficient of Determination Procedures
#' @keywords internal
.n_factors_sescree <- function(eigen_values = NULL, model = "factors") {
nfac <- nFactors::nSeScree(x = eigen_values, cor = TRUE, model = model)$nFactors
.data_frame(
n_Factors = as.numeric(nfac),
Method = c("Scree (SE)", "Scree (R2)"),
Family = "Scree_SE"
)
}
# EGAnet ------------------------
.n_factors_ega <- function(x = NULL,
cor = NULL,
nobs = NULL,
eigen_values = NULL,
type = "FA") {
# Replace with own correlation matrix
junk <- utils::capture.output(suppressWarnings(suppressMessages(
nfac_glasso <- EGAnet::EGA(cor, n = nobs, model = "glasso", plot.EGA = FALSE)$n.dim # nolint
)))
junk <- utils::capture.output(suppressWarnings(suppressMessages(
nfac_TMFG <- .safe(EGAnet::EGA(cor, n = nobs, model = "TMFG", plot.EGA = FALSE)$n.dim, NA) # nolint
)))
.data_frame(
n_Factors = as.numeric(c(nfac_glasso, nfac_TMFG)),
Method = c("EGA (glasso)", "EGA (TMFG)"),
Family = "EGA"
)
}
# psych ------------------------
#' @keywords internal
.n_factors_parallel <- function(x = NULL,
cor = NULL,
nobs = NULL,
type = "FA") {
# Altnerative version of parralel analysis
# Not used because already included in nFactors
if (tolower(type) %in% c("fa", "factor", "efa")) {
fa <- "fa"
} else {
fa <- "pc"
}
insight::check_if_installed("psych")
out <- psych::fa.parallel(cor, n.obs = nobs, fa = fa, plot = FALSE, fm = "ml")
.data_frame(
n_Factors = as.numeric(stats::na.omit(c(out$nfact, out$ncomp))),
Method = "Parallel",
Family = "psych"
)
}
#' @keywords internal
.n_factors_vss <- function(x = NULL,
cor = NULL,
nobs = NULL,
type = "FA",
rotation = "varimax",
algorithm = "default") {
if (algorithm == "default") {
if (tolower(type) %in% c("fa", "factor", "efa")) {
algorithm <- "minres"
} else {
algorithm <- "pc"
}
}
insight::check_if_installed("psych")
# Compute VSS
vss <- psych::VSS(
cor,
n = ncol(x) - 1,
n.obs = nobs,
rotate = rotation,
fm = algorithm,
plot = FALSE
)
# Format results
stats <- vss$vss.stats
stats$map <- vss$map
stats$n_Factors <- seq_len(nrow(stats))
names(stats) <- gsub("cfit.", "VSS_Complexity_", names(stats))
# Indices
vss_1 <- which.max(stats$VSS_Complexity_1)
vss_2 <- which.max(stats$VSS_Complexity_2)
velicer_MAP <- which.min(stats$map)
BIC_reg <- which.min(stats$BIC)
BIC_adj <- which.min(stats$SABIC)
BIC_reg <- ifelse(length(BIC_reg) == 0, NA, BIC_reg)
BIC_adj <- ifelse(length(BIC_adj) == 0, NA, BIC_adj)
.data_frame(
n_Factors = as.numeric(c(vss_1, vss_2, velicer_MAP, BIC_reg, BIC_adj)),
Method = c("VSS complexity 1", "VSS complexity 2", "Velicer's MAP", "BIC", "BIC (adjusted)"),
Family = c("VSS", "VSS", "Velicers_MAP", "BIC", "BIC")
)
}
#' @keywords internal
.n_factors_fit <- function(x = NULL,
cor = NULL,
nobs = NULL,
type = "FA",
rotation = "varimax",
algorithm = "default",
threshold = 0.1) {
if (algorithm == "default") {
if (tolower(type) %in% c("fa", "factor", "efa")) {
algorithm <- "minres"
} else {
algorithm <- "pc"
}
}
insight::check_if_installed("psych")
rez <- data.frame()
for (n in 1:(ncol(cor) - 1)) {
if (tolower(type) %in% c("fa", "factor", "efa")) {
factors <- tryCatch(
suppressWarnings(
psych::fa(
cor,
nfactors = n,
n.obs = nobs,
rotate = rotation,
fm = algorithm
)
),
error = function(e) NA
)
} else {
factors <- tryCatch(
suppressWarnings(
psych::pca(
cor,
nfactors = n,
n.obs = nobs,
rotate = rotation
)
),
error = function(e) NA
)
}
if (all(is.na(factors))) {
next
}
rmsea <- ifelse(is.null(factors$RMSEA), NA, factors$RMSEA[1])
rmsr <- ifelse(is.null(factors$rms), NA, factors$rms)
crms <- ifelse(is.null(factors$crms), NA, factors$crms)
bic <- ifelse(is.null(factors$BIC), NA, factors$BIC)
tli <- ifelse(is.null(factors$TLI), NA, factors$TLI)
rez <- rbind(
rez,
.data_frame(
n = n,
Fit = factors$fit.off,
TLI = tli,
RMSEA = rmsea,
RMSR = rmsr,
CRMS = crms,
BIC = bic
)
)
}
# For fit indices that constantly increase / decrease, we need to find
# an "elbow"/"knee". Here we take the first value that reaches 90 percent
# of the range between the max and the min (when 'threshold = 0.1').
# Fit
if (all(is.na(rez$Fit))) {
fit_off <- NA
} else {
target <- max(rez$Fit, na.rm = TRUE) - threshold * diff(range(rez$Fit, na.rm = TRUE))
fit_off <- rez[!is.na(rez$Fit) & rez$Fit >= target, "n"][1]
}
# TLI
if (all(is.na(rez$TLI))) {
TLI <- NA
} else {
target <- max(rez$TLI, na.rm = TRUE) - threshold * diff(range(rez$TLI, na.rm = TRUE))
TLI <- rez[!is.na(rez$TLI) & rez$TLI >= target, "n"][1]
}
# RMSEA
if (all(is.na(rez$RMSEA))) {
RMSEA <- NA
} else {
target <- min(rez$RMSEA, na.rm = TRUE) + threshold * diff(range(rez$RMSEA, na.rm = TRUE))
RMSEA <- rez[!is.na(rez$RMSEA) & rez$RMSEA <= target, "n"][1]
}
# RMSR
if (all(is.na(rez$RMSR))) {
RMSR <- NA
} else {
target <- min(rez$RMSR, na.rm = TRUE) + threshold * diff(range(rez$RMSR, na.rm = TRUE))
RMSR <- rez[!is.na(rez$RMSR) & rez$RMSR <= target, "n"][1]
}
# CRMS
if (all(is.na(rez$CRMS))) {
CRMS <- NA
} else {
target <- min(rez$CRMS, na.rm = TRUE) + threshold * diff(range(rez$CRMS, na.rm = TRUE))
CRMS <- rez[!is.na(rez$CRMS) & rez$CRMS <= target, "n"][1]
}
# BIC (this is a penalized method so we can just take the one that minimizes it)
BayIC <- ifelse(all(is.na(rez$BIC)), NA, rez[!is.na(rez$BIC) & rez$BIC == min(rez$BIC, na.rm = TRUE), "n"])
.data_frame(
n_Factors = c(fit_off, TLI, RMSEA, RMSR, CRMS, BayIC),
Method = c("Fit_off", "TLI", "RMSEA", "RMSR", "CRMS", "BIC"),
Family = c("Fit", "Fit", "Fit", "Fit", "Fit", "Fit")
)
}
# PCDimension ------------------------
#' @keywords internal
.n_factors_PCDimension <- function(x = NULL, type = "PCA") {
# This package is a strict dependency of PCDimension so if users have the
# former they should have it
insight::check_if_installed(c("ClassDiscovery", "PCDimension"))
# Only applies to PCA with full data
if (tolower(type) %in% c("fa", "factor", "efa") || !is.data.frame(x)) {
return(data.frame())
}
# Randomization-Based Methods
rez_rnd <- PCDimension::rndLambdaF(x)
# Broken-Stick
spca <- ClassDiscovery::SamplePCA(t(x))
lambda <- spca@variances[1:(ncol(x) - 1)]
rez_bokenstick <- PCDimension::bsDimension(lambda)
# Auer-Gervini
ag <- PCDimension::AuerGervini(spca)
agfuns <- list(
twice = PCDimension::agDimTwiceMean,
specc = PCDimension::agDimSpectral,
km = PCDimension::agDimKmeans,
km3 = PCDimension::agDimKmeans3,
# tt=PCDimension::agDimTtest, # known to overestimate
# cpm=PCDimension::makeAgCpmFun("Exponential"), # known to overestimate
tt2 = PCDimension::agDimTtest2,
cpt = PCDimension::agDimCPT
)
rez_ag <- PCDimension::compareAgDimMethods(ag, agfuns)
.data_frame(
n_Factors = as.numeric(c(rez_rnd, rez_bokenstick, rez_ag)),
Method = c(
"Random (lambda)", "Random (F)", "Broken-Stick", "Auer-Gervini (twice)",
"Auer-Gervini (spectral)", "Auer-Gervini (kmeans-2)", "AuerGervini (kmeans-3)",
"Auer-Gervini (T)", "AuerGervini (CPT)"
),
Family = "PCDimension"
)
}
# Re-implementation of nBentler in nFactors ------------------------
#' @keywords internal
.nBentler <- function(x,
N,
model = model,
log = TRUE,
alpha = 0.05,
cor = TRUE,
details = TRUE,
...) {
insight::check_if_installed("nFactors")
lambda <- nFactors::eigenComputes(x, cor = cor, model = model, ...)
if (any(lambda < 0)) {
insight::format_error(
"These indices are only valid with a principal component solution. So, only positive eigenvalues are permitted."
)
}
minPar <- c(min(lambda) - abs(min(lambda)) + 0.001, 0.001)
maxPar <- c(max(lambda), stats::lm(lambda ~ I(rev(seq_along(lambda))))$coef[2])
n <- N
significance <- alpha
min.k <- 3
LRT <- .data_frame(
q = numeric(length(lambda) - min.k), k = numeric(length(lambda) - min.k),
LRT = numeric(length(lambda) - min.k), a = numeric(length(lambda) - min.k),
b = numeric(length(lambda) - min.k),
p = numeric(length(lambda) - min.k),
convergence = numeric(length(lambda) - min.k)
)
bentler.n <- 0
for (i in 1:(length(lambda) - min.k)) {
temp <-
nFactors::bentlerParameters(
x = lambda,
N = n,
nFactors = i,
log = log,
cor = cor,
minPar = minPar,
maxPar = maxPar,
graphic = FALSE
)
LRT[i, 3] <- temp$lrt
LRT[i, 4] <- ifelse(is.null(temp$coef[1]), NA, temp$coef[1])
LRT[i, 5] <- ifelse(is.null(temp$coef[2]), NA, temp$coef[2])
LRT[i, 6] <- ifelse(is.null(temp$p.value), NA, temp$p.value)
LRT[i, 7] <- ifelse(is.null(temp$convergence), NA, temp$convergence)
LRT[i, 2] <- i
LRT[i, 1] <- length(lambda) - i
}
# LRT <- LRT[order(LRT[,1],decreasing = TRUE),]
for (i in 1:(length(lambda) - min.k)) {
if (i == 1) bentler.n <- bentler.n + as.numeric(LRT$p[i] <= significance)
if (i > 1 && LRT$p[i - 1] <= 0.05) {
bentler.n <- bentler.n + as.numeric(LRT$p[i] <= significance)
}
}
if (bentler.n == 0) {
bentler.n <- length(lambda)
}
if (isTRUE(details)) {
details <- LRT
} else {
details <- NULL
}
res <- list(detail = details, nFactors = bentler.n)
class(res) <- c("nFactors", "list")
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.