library(psych)
library(GPArotation)

data(bfi, package = "psych")
X <- na.omit(bfi[, 1:25])        # 25 items; complete cases only

Xs <- scale(X, center = TRUE, scale = TRUE)   # z-standardise
R  <- cor(Xs)                                 # Pearson correlation matrix

# Kaiser-Meyer-Olkin measure of sampling adequacy (overall and per variable)
kmo <- KMO(R)
kmo$MSA          # overall MSA
kmo$MSAi         # per-variable MSA

# Bartlett's test of sphericity
bart <- cortest.bartlett(R, n = nrow(Xs))
bart$chisq; bart$df; bart$p.value

# Kaiser criterion: eigenvalues of the correlation matrix greater than 1
ev <- eigen(R)$values
sum(ev > 1)

# Parallel analysis (Horn, 1965)
pa <- fa.parallel(R, n.obs = nrow(Xs), fa = "both", plot = FALSE)
pa$nfact

# Very Simple Structure and Velicer's Minimum Average Partial
vss <- VSS(Xs, n = 8, rotate = "varimax", fm = "pa", plot = FALSE, SMC = FALSE)
which.min(vss$map)          # factor count minimising MAP
vss$vss.stats$cfit.1        # VSS complexity 1
vss$vss.stats$cfit.2        # VSS complexity 2

nf  <- 5                                        # retained factors
fit <- fa(r        = Xs,
          nfactors = nf,
          fm       = "pa",                      # Principal Axis extraction
          rotate   = "varimax",                 # orthogonal rotation
          scores   = "Bartlett",
          SMC      = FALSE)

L <- as.matrix(fit$loadings)     # rotated factor loadings
fit$communality                  # h2, variance of each variable explained
fit$uniquenesses                 # u2 = 1 - h2

# Variance explained per factor, and cumulatively
prop_var <- colSums(L^2) / nrow(L)
cum_var  <- cumsum(prop_var)

# Each variable is assigned to the factor carrying its largest loading,
# provided that loading is at least 0.4
dominant <- apply(abs(L), 1, function(z)
  if (max(z) >= 0.4) which.max(z) else NA_integer_)

alphas <- sapply(seq_len(nf), function(f) {
  items <- rownames(L)[which(dominant == f)]
  if (length(items) < 2) return(NA_real_)      # alpha needs >= 2 items
  psych::alpha(X[, items, drop = FALSE], check.keys = TRUE)$total$raw_alpha
})
alphas

fit_ob <- fa(r = Xs, nfactors = nf, fm = "pa",
             rotate = "oblimin", scores = "Bartlett", SMC = FALSE)

fit_ob$Phi        # factor correlation matrix; NULL under orthogonal rotation

fit$dof                 # degrees of freedom
fit$rms                 # root mean square residual (RMSR)
fit$RMSEA               # RMSEA and its confidence bounds
fit$TLI                 # Tucker-Lewis index
fit$STATISTIC; fit$PVAL # chi-square and its p-value
fit$BIC

fit$model               # reproduced (model-implied) correlation matrix
fit$residual            # observed minus reproduced

scores <- fit$scores    # Bartlett factor scores, one row per observation
head(scores)