lp.RdPerforms \(L^p\) rotation to obtain sparse loadings.
Initial factor loadings matrix to be rotated.
Initial rotation matrix.
Component-wise \(L^p\) where 0 < p \(=<\) 1.
Not recommended for \(L^p\) rotation.
Convergence is assumed when the norm of the gradient is smaller than eps.
Maximum number of iterations allowed in the main loop.
Maximum iterations for GPA rotation. The goal is to decrease the objective value, not optimize the inner loop. Warnings may appear, but they can be ignored if the main loop converges.
A GPArotation object, which is a list containing:
Rotated loadings matrix, with one column per factor. If randomStarts
were used, this contains the loadings with the lowest criterion value.
Rotation matrix, satisfying loadings %*% t(Th) = A.
Matrix recording iteration details during optimization.
String indicating the rotation objective function.
Logical indicating whether the rotation is orthogonal.
Logical indicating whether convergence was achieved. Convergence is controlled element-wise by tolerance.
Covariance matrix of rotated factors, t(Th) %*% Th.
These functions optimize an \(L^p\) rotation objective, where 0 < p =< 1. A smaller p promotes sparsity in the loading matrix but increases computational difficulty. For guidance on choosing p, see the Concluding Remarks in the references.
Since the \(L^p\) function is nonsmooth, a different optimization method is required compared to smooth rotation criteria. Two new functions, GPForth.lp and GPFoblq.lp, replace GPForth and GPFoblq for orthogonal and oblique \(L^p\) rotations, respectively.
The optimization method follows an iterative reweighted least squares (IRLS) approach. It approximates the nonsmooth objective function with a smooth weighted least squares function in the main loop and optimizes it using GPA in the inner loop.
Normalization is not recommended for \(L^p\) rotation. Its use may have unexpected results.
lpT,
lpQ,
vgQ.lp.wls
Liu, X., Wallin, G., Chen, Y., & Moustaki, I. (2023). Rotation to sparse loadings using \(L^p\) losses and related inference problems. Psychometrika, 88(2), 527–553.
data("WansbeekMeijer", package = "GPArotation")
fa.unrotated <- factanal(factors = 2, covmat = NetherlandsTV, rotation = "none")
options(warn = -1)
# Orthogonal rotation
# Single start from random position
fa.lpT1 <- GPForth.lp(loadings(fa.unrotated), p = 1)
# 10 random starts
fa.lpT <- lpT(loadings(fa.unrotated), Tmat=Random.Start(2), p = 1, randomStarts = 10)
print(fa.lpT, digits = 5, sortLoadings = FALSE, Table = TRUE, rotateMat = TRUE)
#> Orthogonal rotation method Lp rotation, p=1 converged at lowest minimum.
#> Of 10 random starts 100% converged, 50% at the same lowest minimum.
#> Random starts converged to 2 different local minima
#> Loadings at lowest minimum:
#> Factor1 Factor2
#> NL1 0.66539031 -0.42787
#> TV2 0.65607038 -0.52479
#> NL3 0.64841375 -0.42123
#> RTL4 0.11492304 -0.69264
#> RTL5 0.07431037 -0.75643
#> Veronica -0.00980093 -0.81783
#> SBS6 0.00046327 -0.72669
#>
#> Factor1 Factor2
#> SS loadings 1.31244 2.88478
#> Proportion Var 0.18749 0.41211
#> Cumulative Var 0.18749 0.59960
#>
#> Rotating matrix:
#> [,1] [,2]
#> [1,] 0.48590 -0.87401
#> [2,] -0.87401 -0.48590
#>
#> Iteration table:
#> iter f time
#> elapsed 20 0.93384 0.017
p <- 1
# Oblique rotation
# Single start
fa.lpQ1 <- GPFoblq.lp(loadings(fa.unrotated), p = p)
# 10 random starts
fa.lpQ <- lpQ(loadings(fa.unrotated), p = p, randomStarts = 10)
summary(fa.lpQ, Structure = TRUE)
#> Oblique rotation method Lp rotation, p=1 converged in 13 iterations.
#> Pattern (loadings):
#> Factor1 Factor2
#> NL1 0.000863 -0.79155
#> TV2 -0.101454 -0.78160
#> NL3 -0.003413 -0.77140
#> RTL4 -0.614923 -0.14354
#> RTL5 -0.704348 -0.09622
#> Veronica -0.819252 0.00257
#> SBS6 -0.722079 -0.00857
#>
#> Structure:
#> Factor1 Factor2
#> NL1 -0.422 -0.791
#> TV2 -0.519 -0.836
#> NL3 -0.415 -0.773
#> RTL4 -0.692 -0.472
#> RTL5 -0.756 -0.472
#> Veronica -0.818 -0.435
#> SBS6 -0.727 -0.394
# this functions ensures consistent ordering of factors of a
# GPArotation object for cleaner comparison
# Inspired by fungible::orderFactors and fungible::faSort functions
sortFac <- function(x){
# Only works for object of class GPArotation
if (!inherits(x, "GPArotation")) {stop("argument not of class GPArotation")}
cln <- colnames(x$loadings)
# ordering for oblique slightly different from orthogonal
ifelse(x$orthogonal, vx <- colSums(x$loadings^2),
vx <- diag(x$Phi %*% t(x$loadings) %*% x$loadings) )
# sort by squared loadings from high to low
vxo <- order(vx, decreasing = TRUE)
vx <- vx[vxo]
# maintain the right sign
Dsgn <- diag(sign(colSums(x$loadings^3))) [ , vxo]
x$Th <- x$Th %*% Dsgn
x$loadings <- x$loadings %*% Dsgn
if (match("Phi", names(x))) {
# If factor is negative, reverse corresponding factor correlations
x$Phi <- t(Dsgn) %*% x$Phi %*% Dsgn
}
colnames(x$loadings) <- cln
x
}
# seed set to see the results of sorting
set.seed(1020)
fa.lpQ1 <- lpQ(loadings(fa.unrotated),p=1,randomStarts=10)
fa.lpQ0.5 <- lpQ(loadings(fa.unrotated),p=0.5,randomStarts=10)
fa.geo <- geominQ(loadings(fa.unrotated), randomStarts=10)
# with ordered factor loadings
res <- round(cbind(sortFac(fa.lpQ1)$loadings, sortFac(fa.lpQ0.5)$loadings,
sortFac(fa.geo)$loadings),3)
print(c("oblique- Lp 1 Lp 0.5 geomin")); print(res)
#> [1] "oblique- Lp 1 Lp 0.5 geomin"
#> Factor1 Factor2 Factor1 Factor2 Factor1 Factor2
#> NL1 -0.001 0.792 -0.002 0.792 -0.020 0.802
#> TV2 0.101 0.782 0.101 0.781 0.085 0.788
#> NL3 0.003 0.771 0.002 0.772 -0.015 0.782
#> RTL4 0.615 0.144 0.618 0.138 0.627 0.119
#> RTL5 0.704 0.096 0.708 0.089 0.719 0.067
#> Veronica 0.819 -0.003 0.824 -0.011 0.839 -0.038
#> SBS6 0.722 0.009 0.726 0.002 0.740 -0.023
# without ordered factor loadings
res <- round(cbind(fa.lpQ1$loadings, fa.lpQ0.5$loadings, fa.geo$loadings),3)
print(c("oblique- Lp 1 Lp 0.5 geomin")); print(res)
#> [1] "oblique- Lp 1 Lp 0.5 geomin"
#> Factor1 Factor2 Factor1 Factor2 Factor1 Factor2
#> NL1 0.792 -0.001 -0.792 0.002 0.020 -0.802
#> TV2 0.782 0.101 -0.781 -0.101 -0.085 -0.788
#> NL3 0.771 0.003 -0.772 -0.002 0.015 -0.782
#> RTL4 0.144 0.615 -0.138 -0.618 -0.627 -0.119
#> RTL5 0.096 0.704 -0.089 -0.708 -0.719 -0.067
#> Veronica -0.003 0.819 0.011 -0.824 -0.839 0.038
#> SBS6 0.009 0.722 -0.002 -0.726 -0.740 0.023
options(warn = 0)