Comment calculer les composants principaux à rotation varimax dans R?

J'ai exécuté PCA sur 25 variables et sélectionné les 7 meilleurs PC utilisant prcomp.

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

J'ai ensuite fait une rotation varimax sur ces composants.

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

Et maintenant, je souhaite faire pivoter varimax les données pivotées par PCA (car elles ne font pas partie de l'objet varimax - seulement la matrice des chargements et la matrice de rotation). J'ai lu que pour ce faire on multiplie la transposition de la matrice de rotation par la transposition des données donc j'aurais fait ça:

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

Mais cela n'a pas de sens car les dimensions de la matrice transposée ci-dessus sont respectivement $7\times 7$ et $7 \times 16933$ et il me restera donc une matrice de seulement $7$ lignes, plutôt que $16933$ lignes ... est-ce que quelqu'un sait ce que je fais mal ici ou quelle devrait être ma dernière ligne? Dois-je simplement transposer en arrière après?

"Rotations" est une approche développée en analyse factorielle; là des rotations (comme par exemple varimax) sont appliquées aux chargements , pas aux vecteurs propres de la matrice de covariance. Les charges sont des vecteurs propres mis à l'échelle par les racines carrées des valeurs propres respectives. Après la rotation varimax, les vecteurs de chargement ne sont plus orthogonaux (même si la rotation est appelée "orthogonale"), donc on ne peut pas simplement calculer des projections orthogonales des données sur les directions de chargement tournées.

@ La réponse de FTusell suppose que la rotation varimax est appliquée aux vecteurs propres (pas aux chargements). Ce serait assez peu conventionnel. Veuillez consulter mon compte rendu détaillé de PCA + varimax pour plus de détails: PCA suivi d'une rotation (comme varimax) est-il toujours PCA? En bref, si nous regardons la SVD de la matrice de données $X=USV^\top$ , alors faire tourner les chargements signifie insérer $RR^\top$ pour une matrice de rotation $R$ comme suit: $X=(UR)(R^\top SV^\top).$

Si la rotation est appliquée aux chargements (comme c'est généralement le cas), il existe au moins trois façons simples de calculer les PC à rotation varimax dans R:

1. Ils sont facilement disponibles via la fonction psych::principal(démontrant qu'il s'agit bien de l'approche standard). Notez qu'il renvoie des scores standardisés , c'est-à-dire que tous les PC ont une variance d'unité.

2. On peut utiliser manuellement la varimaxfonction pour faire pivoter les chargements, puis utiliser les nouveaux chargements tournés pour obtenir les scores; il faut multiplier les données avec le pseudo-inverse transposé des chargements tournés (voir formules dans cette réponse de @ttnphns ). Cela produira également des scores standardisés.

3. On peut utiliser la varimaxfonction pour faire pivoter les chargements, puis utiliser la $rotmatmatrice de rotation pour faire pivoter les scores standardisés obtenus avec prcomp.

Les trois méthodes donnent le même résultat:

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

Cela donne trois sorties identiques:

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

---

Remarque: La varimaxfonction dans R utilise des normalize = TRUE, eps = 1e-5paramètres par défaut ( voir la documentation ). On pourrait vouloir changer ces paramètres (diminuer la epstolérance et prendre soin de la normalisation de Kaiser) lors de la comparaison des résultats avec d'autres logiciels tels que SPSS. Je remercie @GottfriedHelms d'avoir porté cela à mon attention. [Remarque: ces paramètres fonctionnent lorsqu'ils sont transmis à la varimaxfonction, mais ne fonctionnent pas lorsqu'ils sont transmis à la psych::principalfonction. Cela semble être un bug qui sera corrigé.]

Vous devez utiliser la matrice $loadings, pas $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 matrice $rotmatest la matrice orthogonale qui produit les nouveaux chargements à partir de ceux non tournés.

EDIT au 12 février 2015:

Comme indiqué à juste titre ci-dessous par @amoeba (voir également son article précédent ainsi qu'un autre article de @ttnphns ), cette réponse n'est pas correcte. Considérons un matrice de données X . La décomposition en valeurs singulières est X = U S V T où V a pour colonnes (normalisé) les vecteurs propres de X ' X . Maintenant, une rotation est un changement de coordonnées et revient à écrire l'égalité ci-dessus comme: X = ( U S T ) ( T T V T$n\times m$$X$

$$X = USV^T$$ $V$$X'X$ ≈ avec T étant une matrice orthogonale choisie pour obtenir un V ∗ proche de clairsemé (contraste maximum entre les entrées, en gros). Maintenant,si c'était tout, ce qui n'est pas le cas, on pourrait post-multiplier l'égalité ci-dessus par V ∗ pour obtenir les scores U ∗ comme X ( V ∗ ) T , Mais bien sûr, nous ne faisons jamais tourner tous les PC. Nous considérons plutôt un sous-ensemble de k < m qui fournit toujours uneapproximationdécente de rang k de X , X ( $$X = (UST)(T^TV^T) = U^*V^*$$ $T$$V^*$$V^*$$U^*$$X(V^*)^T$$k<m$$k$$X$ donc la solution tournée est maintenant X ≈ ( U k S k T k ) ( T T k V T k ) = U ∗ k V ∗ k où maintenant V ∗ k est un k × n matrice. On ne peut plus simplement multiplier X par la transposition 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^*$, mais nous devons plutôt recourir à l'une des solutions décrites par @amoeba.

En d'autres termes, la solution que j'ai proposée n'est correcte que dans le cas particulier où elle serait inutile et absurde.

Un grand merci à @amoeba pour m'avoir clarifié cette question; Je vis avec cette idée fausse depuis des années.

$S$$V$$L$$V$$S$$v_i^TX$ $(i=1,\ldots,m)$$\|v_i\|=1$. De toute façon est acceptable, je pense, et tout le reste (comme dans l'analyse biplot).

NOUVELLE MODIFICATION 12 février 2015

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

Je cherchais une solution qui fonctionne pour PCA réalisée en utilisant ade4 .

Veuillez trouver la fonction ci-dessous:

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

^{Créé le 2020-01-14 par le package reprex (v0.3.0)}

J'espère que cette aide!