¿Cómo calcular componentes principales rotados con varimax en R?

Ejecuté PCA en 25 variables y seleccioné las 7 mejores PC usando prcomp.

prc <- prcomp(pollutions, center=T, scale=T, retx=T)

Luego hice rotación varimax en esos componentes.

varimax7 <- varimax(prc$rotation[,1:7])

Y ahora deseo varimax rotar los datos rotados por PCA (ya que no es parte del objeto varimax, solo la matriz de cargas y la matriz de rotación). Leí que para hacer esto, multiplique la transposición de la matriz de rotación por la transposición de los datos, así que habría hecho esto:

newData <- t(varimax7$rotmat) %*% t(prc$x[,1:7])

Pero eso no tiene sentido ya que las dimensiones de las transposiciones de la matriz anteriores son $7\times 7$ y $7 \times 16933$ respectivamente, por lo que me quedará con una matriz de solo $7$ filas, en lugar de filas ... ¿Alguien sabe lo que yo estoy haciendo mal aquí o cuál debería ser mi línea final? ¿Solo necesito transponerme después?$16933$

"Rotaciones" es un enfoque desarrollado en el análisis factorial; allí las rotaciones (como, por ejemplo, varimax) se aplican a las cargas , no a los vectores propios de la matriz de covarianza. Las cargas son vectores propios escalados por las raíces cuadradas de los valores propios respectivos. Después de la rotación varimax, los vectores de carga ya no son ortogonales (aunque la rotación se llama "ortogonal"), por lo que no se puede calcular simplemente las proyecciones ortogonales de los datos en las direcciones de carga rotadas.

La respuesta de @ FTusell supone que la rotación varimax se aplica a los vectores propios (no a las cargas). Esto sería bastante poco convencional. Consulte mi cuenta detallada de PCA + varimax para obtener más detalles: ¿PCA seguido de una rotación (como varimax) sigue siendo PCA? Brevemente, si miramos el SVD de la matriz de datos , entonces rotar las cargas significa insertar R R ⊤ para alguna matriz de rotación R de la siguiente manera: X = ( U R ) ( R ⊤ S V ⊤ ) $X=USV^\top$$RR^\top$$R$$X=(UR)(R^\top SV^\top).$

Si la rotación se aplica a las cargas (como suele serlo), entonces hay al menos tres formas fáciles de calcular las PC rotadas con varimax en R:

1. Están fácilmente disponibles a través de la función psych::principal(lo que demuestra que este es realmente el enfoque estándar). Tenga en cuenta que devuelve puntajes estandarizados , es decir, todas las PC tienen varianza unitaria.

2. Se puede usar manualmente la varimaxfunción para rotar las cargas, y luego usar las nuevas cargas rotadas para obtener los puntajes; uno necesita multiplicar los datos con el pseudoinverso inverso de las cargas rotadas (vea las fórmulas en esta respuesta de @ttnphns ). Esto también generará puntajes estandarizados.

3. Se puede usar la varimaxfunción para rotar las cargas, y luego usar la $rotmatmatriz de rotación para rotar los puntajes estandarizados obtenidos con prcomp.

Los tres métodos producen el mismo resultado:

irisX <- iris[,1:4]      # Iris data
ncomp <- 2

pca_iris_rotated <- psych::principal(irisX, rotate="varimax", nfactors=ncomp, scores=TRUE)
print(pca_iris_rotated$scores[1:5,])  # Scores returned by principal()

pca_iris        <- prcomp(irisX, center=T, scale=T)
rawLoadings     <- pca_iris$rotation[,1:ncomp] %*% diag(pca_iris$sdev, ncomp, ncomp)
rotatedLoadings <- varimax(rawLoadings)$loadings
invLoadings     <- t(pracma::pinv(rotatedLoadings))
scores          <- scale(irisX) %*% invLoadings
print(scores[1:5,])                   # Scores computed via rotated loadings

scores <- scale(pca_iris$x[,1:2]) %*% varimax(rawLoadings)$rotmat
print(scores[1:5,])                   # Scores computed via rotating the scores

Esto produce tres salidas idénticas:

1 -1.083475  0.9067262
2 -1.377536 -0.2648876
3 -1.419832  0.1165198
4 -1.471607 -0.1474634
5 -1.095296  1.0949536

---

Nota: La varimaxfunción en R usa normalize = TRUE, eps = 1e-5parámetros por defecto ( ver documentación ). Es posible que desee cambiar estos parámetros (disminuir la epstolerancia y ocuparse de la normalización de Kaiser) al comparar los resultados con otro software como SPSS. Agradezco a @GottfriedHelms por llamar mi atención sobre esto. [Nota: estos parámetros funcionan cuando se pasan a la varimaxfunción, pero no funcionan cuando se pasan a la psych::principalfunción. Esto parece ser un error que se solucionará.]

Necesita usar la matriz $loadings, no $rotmat:

 x <- matrix(rnorm(600),60,10)
 prc <- prcomp(x, center=TRUE, scale=TRUE)
 varimax7 <- varimax(prc$rotation[,1:7])
 newData <- scale(x) %*% varimax7$loadings

La matriz $rotmates la matriz ortogonal que produce las nuevas cargas de los no rotados.

EDITAR a partir del 12 de febrero de 2015:

Como señala acertadamente a continuación @amoeba (consulte también su publicación anterior y otra publicación de @ttnphns ), esta respuesta no es correcta. Considere un matriz de datos X . La descomposición de valor singular es X = U S V T donde V tiene como sus columnas los (normalizado) vectores propios de X ' X . Ahora, una rotación es un cambio de coordenadas y equivale a escribir la igualdad anterior como: X = ( U S T ) ( T T V T$n\times m$$X$

$$X = USV^T$$ $V$$X'X$ con T siendo una matriz ortogonal elegida para lograr un V ∗ cercano a disperso (contraste máximo entre entradas, hablando en términos generales). Ahora,si eso fuera todo, lo que no es, uno podría multiplicar posteriormente la igualdad anterior por V ∗ para obtener puntuaciones U ∗ como X ( V ∗ ) T , pero por supuesto, nunca rotamos todas las PC. Por el contrario, consideramos un subconjunto de k < m que proporciona todavía unaaproximacióndecente de rango k de X , X $$X = (UST)(T^TV^T) = U^*V^*$$ $T$$V^*$$V^*$$U^*$$X(V^*)^T$$k<m$$k$$X$ entonces la solución rotada es ahora X ≈ ( U k S k T k ) ( T T k V T k ) = U ∗ k V ∗ k donde ahora V ∗ k es una k × n matriz. Ya no podemos simplemente multiplicar X por la transposición de V ∗ $$X \approx (U_kS_k)(V_k^T)$$ $$X \approx (U_kS_kT_k)(T_k^TV_k^T) = U_k^*V_k^*$$ $V_k^*$$k\times n$$X$$V_k^*$, sino que necesitamos recurrir a una de las soluciones descritas por @amoeba.

En otras palabras, la solución que propuse solo es correcta en el caso particular en el que sería inútil y sin sentido.

Agradezco sinceramente a @amoeba por aclararme este asunto; He estado viviendo con este concepto erróneo durante años.

$S$$V$$L$$V$$S$$v_i^TX$ $(i=1,\ldots,m)$$\|v_i\|=1$. De cualquier manera es aceptable, creo, y todo lo demás (como en el análisis biplot).

MÁS EDITAR 12 de febrero de 2015

$V_k^*$$V_k^*(V_k^*)^T$$X(V_k^*)^T \approx U_k^*$

Estaba buscando una solución que funcione para PCA realizada con ade4 .

Encuentre la función a continuación:

library(ade4)

irisX <- iris[,1:4]      # Iris data
ncomp <- 2
# With ade4
dudi_iris <- dudi.pca(irisX, scannf = FALSE, nf = ncomp)

rotate_dudi.pca <- function(pca, ncomp = 2) {

  rawLoadings <- as.matrix(pca$c1[,1:ncomp]) %*% diag(sqrt(pca$eig), ncomp, ncomp)
  pca$c1 <- rawLoadings
  pca$li <- scale(pca$li[,1:ncomp]) %*% varimax(rawLoadings)$rotmat

  return(pca)
} 
rot_iris <- rotate_dudi.pca(pca = dudi_iris, ncomp = ncomp)
print(rot_iris$li[1:5,])                   # Scores computed via rotating the scores
#>        [,1]       [,2]
#> 1 -1.083475 -0.9067262
#> 2 -1.377536  0.2648876
#> 3 -1.419832 -0.1165198
#> 4 -1.471607  0.1474634
#> 5 -1.095296 -1.0949536

^{Creado el 14-01-2020 por el paquete reprex (v0.3.0)}

¡Espero que esto ayude!