GPArotation (version 2025.3-1)

lp: \(L^p\) Rotation

Description

Performs \(L^p\) rotation to obtain sparse loadings.

Usage

GPForth.lp(A, Tmat=diag(rep(1, ncol(A))), p=1, normalize=FALSE, eps=1e-05, 
            maxit=10000, gpaiter=5) 
     GPFoblq.lp(A, Tmat=diag(rep(1, ncol(A))), p=1, normalize=FALSE, eps=1e-05, 
            maxit=10000, gpaiter=5)

Value

A GPArotation object, which is a list containing:

loadings

Rotated loadings matrix, with one column per factor. If randomStarts were used, this contains the loadings with the lowest criterion value.

Th

Rotation matrix, satisfying loadings %*% t(Th) = A.

Table

Matrix recording iteration details during optimization.

method

String indicating the rotation objective function.

orthogonal

Logical indicating whether the rotation is orthogonal.

convergence

Logical indicating whether convergence was achieved. Convergence is controlled element-wise by tolerance.

Phi

Covariance matrix of rotated factors, t(Th) %*% Th.

Arguments

A

Initial factor loadings matrix to be rotated.

Tmat

Initial rotation matrix.

p

Component-wise \(L^p\) where 0 < p \(=<\) 1.

normalize

Not recommended for \(L^p\) rotation.

eps

Convergence is assumed when the norm of the gradient is smaller than eps.

maxit

Maximum number of iterations allowed in the main loop.

gpaiter

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.

Author

Xinyi Liu, with minor modifications for GPArotation by C. Bernaards

Details

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.

References

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.

See Also

lpT, lpQ, vgQ.lp.wls

Examples

  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)
  
  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)


  # 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)

  # 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)

  options(warn = 0)