Wie berechnet man varimaxgedrehte Hauptkomponenten in R?

Ich habe PCA mit 25 Variablen ausgeführt und die 7 besten PCs mit ausgewählt prcomp.

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

Ich habe dann Varimax-Rotation an diesen Komponenten durchgeführt.

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

Und jetzt möchte ich die PCA-gedrehten Daten mit Varimax drehen (da sie nicht Teil des Varimax-Objekts sind - nur die Ladematrix und die Rotationsmatrix). Ich habe gelesen, dass Sie dazu die Transponierung der Rotationsmatrix mit der Transponierung der Daten multiplizieren, also hätte ich dies getan:

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

Das macht aber keinen Sinn, da die Dimensionen der oben transponierten Matrix $7\times 7$ bzw. $7 \times 16933$ sind und ich nur noch eine Matrix mit $7$ Zeilen anstelle von $16933$ Zeilen habe ... weiß jemand, was ich mache ich hier falsch oder was soll meine letzte Zeile sein? Muss ich danach nur noch zurück transponieren?

r pca factor-rotation Scott
quelle

Antworten:

"Rotationen" ist ein Ansatz, der in der Faktorenanalyse entwickelt wurde. dort werden Drehungen (wie z. B. Varimax) auf Ladungen angewendet , nicht auf Eigenvektoren der Kovarianzmatrix. Belastungen sind Eigenvektoren, die durch die Quadratwurzeln der jeweiligen Eigenwerte skaliert werden. Nach der Varimax-Drehung sind die Ladevektoren nicht mehr orthogonal (obwohl die Drehung "orthogonal" genannt wird), so dass man nicht einfach orthogonale Projektionen der Daten auf die gedrehten Laderichtungen berechnen kann.

Bei der Antwort von @ FTusell wird davon ausgegangen, dass die Varimax-Rotation auf die Eigenvektoren angewendet wird (nicht auf Lasten). Das wäre ziemlich unkonventionell. Weitere Informationen finden Sie in meinem ausführlichen Bericht zu PCA + varimax: Wird nach PCA noch eine Rotation (z. B. varimax) durchgeführt? Betrachten wir kurz die SVD der Datenmatrix , so bedeutet das Drehen der Ladungen das Einfügen von für eine Rotationsmatrix wie folgt: $X=USV^\top$ $RR^\top$ $R$ $X=(UR)(R^\top SV^\top).$

Wenn Ladungen gedreht werden (wie es normalerweise der Fall ist), gibt es mindestens drei einfache Methoden, um varimax-gedrehte PCs in R zu berechnen:

Sie sind über die Funktion leicht verfügbar psych::principal(was zeigt, dass dies tatsächlich der Standardansatz ist). Beachten Sie, dass es standardisierte Werte zurückgibt , dh alle PCs haben eine Einheitsvarianz.
Man kann die varimaxFunktion manuell verwenden , um die Ladungen zu drehen, und dann die neuen gedrehten Ladungen verwenden, um die Bewertungen zu erhalten; man muss die Daten mit der transponierten Pseudoinverse der gedrehten Ladungen multiplizieren (siehe Formeln in dieser Antwort von @ttnphns ). Dies wird auch standardisierte Ergebnisse liefern.
Man kann die varimaxFunktion verwenden, um die Ladungen zu drehen, und dann die $rotmatRotationsmatrix verwenden, um die mit erhaltenen standardisierten Punktzahlen zu drehen prcomp.

Alle drei Methoden führen zu demselben Ergebnis:

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

Dies ergibt drei identische Ausgaben:

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

Hinweis: Die varimaxFunktion in R verwendet normalize = TRUE, eps = 1e-5standardmäßig Parameter ( siehe Dokumentation ). Möglicherweise möchten Sie diese Parameter ändern (verringern Sie die epsToleranz und achten Sie auf die Kaiser-Normalisierung), wenn Sie die Ergebnisse mit anderen Programmen wie SPSS vergleichen. Ich danke @GottfriedHelms, dass er mich darauf aufmerksam gemacht hat. [Hinweis: Diese Parameter funktionieren, wenn sie an die varimaxFunktion übergeben werden, aber nicht, wenn sie an die psych::principalFunktion übergeben werden. Dies scheint ein Fehler zu sein, der behoben wird.]

Amöbe sagt Reinstate Monica
quelle

Ich sehe das jetzt und ich denke, du hast recht. Ich bearbeite meine ursprüngliche Antwort (oder füge eine weitere hinzu), um die Ursache der Diskrepanz zu ermitteln. Ich mochte Ihre und @ttnphns sehr vollständige und aufschlussreiche Antworten, die detaillierte Erklärungen lieferten, die normalerweise nicht in Büchern zu finden sind.

F. Tusell

@amoeba Ich versuche, eine PCA + Varimax mit zu tun principal , prcompund princomp, aber die daraus resultierenden Belastungen / Studie Schlussfolgerungen sind sehr verschieden voneinander. Nach meinem Verständnis geben prcomp und princomp keine standardisierten Scores oder Ladungen zurück. Meine Frage ist: Was ist der beste Ansatz? Will ich wirklich standardisierte Ergebnisse? Ist das nicht mein Code , pca_iris <- prcomp(irisX, center=T, scale=T)gefolgt von varimax(pca_iris$rotation)$loadingsso richtig wie Sie oben?

JMarcelino

@JMarcelino, nein, dein Code dreht Varimax auf den Eigenvektoren, nicht auf den Ladungen. So wird Varimax-Rotation normalerweise nicht verstanden oder angewendet.

Amöbe sagt Reinstate Monica

@JMarcelino, fragst du, warum die Mathematik so funktioniert, wie ich es in Methode 2 gesagt habe? Es ist einfach, wenn Sie mit dieser Art linearer Algebra vertraut sind. PCA ist SVD-Zerlegung

. Eine Rotation wie varimax anzuwenden bedeutet,

X = U S V^{⊤}

$X=USV^\top$

R R^{⊤}

$RR^\top$ für eine Rotationsmatrix

wie folgt

R

$R$

. Gedrehte Lasten sind

X = U R R^{⊤} S V^{⊤}

$X=URR^\top SV^\top$

, rotierte Standardwerte sind

L = V S R / \sqrt{n - 1}

$L=VSR/\sqrt{n-1}$

, also

Du kennst

und

; Wie komme ich nach

? Nun, die Antwort lautet

T = U R \sqrt{n - 1}

$T=UR\sqrt{n-1}$

X = T L^{⊤} .

$X=TL^\top.$

X

$X$

L

$L$

T

$T$

T = X (L^{⊤})^{+} = X (L^{+})^{⊤} .

$T=X(L^\top)^+ = X(L^+)^\top.$

Amöbe sagt Reinstate Monica

Ich erhielt eine Antwort vom Betreuer des Pakets, Prof. Revelle. Es scheint ein Fehler im Umgang mit den Parametern im zu seinprincipal Prozedur zu sein, die immer mit Kaiser-Normalisierung und eps = 1e-5 rechnet. Bisher gibt es keine Informationen, warum auf r-fiddle.org die Version korrekt funktioniert. Wir sollten also auf Updates warten - und alle jetzt veralteten Kommentare löschen. Amöbe - es wäre gut, die Bemerkung in Ihrer Antwort entsprechend zu aktualisieren. Vielen Dank für die Zusammenarbeit!

Gottfried Helms

Sie müssen die Matrix verwenden $loadings, nicht $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

Die Matrix $rotmatist die orthogonale Matrix, die aus den nicht gedrehten Ladungen die neuen Ladungen erzeugt.

BEARBEITEN ab dem 12. Februar 2015:

$n\times m$ $X$

X = U S V^{T}

$X = USV^T$

V

$V$

X^{'} X

$X'X$

X = (U S T) (T^{T} V^{T}) = U^{*} V^{*}

$X = (UST)(T^TV^T) = U^*V^*$

T

$T$

V^{*}

$V^*$

V^{*}

$V^*$

U^{*}

$U^*$

X (V^{*})^{T}

$X(V^*)^T$

k < m

$k<m$

k

$k$

X

$X$

X \approx (U_{k} S_{k}) (V_{k}^{T})

$X \approx (U_kS_k)(V_k^T)$

X \approx (U_{k} S_{k} T_{k}) (T_{k}^{T} V_{k}^{T}) = U_{k}^{*} V_{k}^{*}

$X \approx (U_kS_kT_k)(T_k^TV_k^T) = U_k^*V_k^*$

V_{k}^{*}

$V_k^*$

k \times n

$k\times n$

X

$X$

V_{k}^{*}

$V_k^*$ , sondern wir müssen auf eine der von @amoeba beschriebenen Lösungen zurückgreifen.

Mit anderen Worten, die von mir vorgeschlagene Lösung ist nur dann richtig, wenn sie unbrauchbar und unsinnig wäre.

Ein herzliches Dankeschön geht an @amoeba, dass sie mir diese Angelegenheit klar gemacht haben. Ich lebe seit Jahren mit diesem Missverständnis.

$S$ $V$ $L$ $V$ $S$ $v_i^TX$ $(i=1,\ldots,m)$ $\|v_i\|=1$ . Ich denke, jeder Weg ist akzeptabel und alles dazwischen (wie bei der Biplot-Analyse).

WEITERE ÄNDERUNGEN 12. Februar 2015

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

F. Tusell
quelle

Ah, richtig großartig. Ich war verwirrt, weil die Ladevorgänge für den Prcomp als "Rotation" bezeichnet werden, hätte die Hilfe besser lesen sollen. Bedeutet das, dass ich meine Daten zentrieren und skalieren sollte, bevor ich sie mit meinen varimax $ loadings multipliziere, da ich "center = TRUE, scale = TRUE" in der prcomp-Methode verwende?

Scott

Ja, guter Punkt, mein Fehler. Das Zentrieren spielt keine Rolle, als ob nur die Punkte verschoben würden, aber der Maßstab sollte derselbe sein, der zur Berechnung der Hauptkomponenten verwendet wird, die für die Skalierung nicht invariant sind.

F. Tusell

Ich habe vergessen zu erwähnen, dass Sie sich die Funktion factanal ansehen möchten, wenn Sie dies noch nicht getan haben. Es wird eher eine Faktorenanalyse als eine Hauptkomponente durchgeführt, die Ergebnisse werden jedoch direkt zurückgegeben.

F. Tusell

-1. Ich glaube, dass diese Antwort nicht korrekt ist und habe meine eigene Antwort veröffentlicht, um sie zu demonstrieren. Man kann keine gedrehten Punkte durch orthogonale Projektion auf die gedrehten Ladungen erhalten (weil sie nicht mehr orthogonal sind). Der einfachste Weg, um die richtigen Ergebnisse zu erhalten, ist die Verwendung von psych::principal. [Abgesehen davon habe ich Ihre Antwort bearbeitet, um die Skalierung einzufügen, wie in den obigen Kommentaren erläutert.]

Amöbe sagt Reinstate Monica

Es tut mir leid. ich meinte

V_{k}^{*}

$V_k^*$ ist

k \times n

$k\times n$ . Ich werde es jetzt korrigieren. Und ... ja, jetzt wo ich es mir ansehe,

V

$V$ hat also orthogonale Spalten

(T_{k}^{T} V_{k}^{T}) (V_{k} T_{k})

$(T_k^TV_k^T)(V_kT_k)$ würde uns trotzdem eine Einheitsmatrix besorgen, oder? Wenn ja, habe ich das Originalplakat nicht irregeführt, du hebst eine Last von meiner Seele!

F. Tusell

Ich war auf der Suche nach einer Lösung, die für PCA mit ade4 funktioniert .

Hier finden Sie die Funktion:

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

^{Erstellt am 2020-01-14 von der Paket reprex (v0.3.0)}

Ich hoffe das hilft!

Alain Danet
quelle

Sie müssen diesen Raum für eine Antwort verwenden.

Michael R. Chernick

Es schien mir, dass es gültig ist, der Vollständigkeit halber eine Antwort hinzuzufügen. Beispiel für diese Frage: stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2 . Gerne verschiebe ich mein Angebot bei Bedarf.

Alain Danet

Ich habe es falsch verstanden, weil es sich anhörte, als würden Sie einen Fehler in einer der Antworten korrigieren. Ich sehe, dass es eine Ergänzung für ein bestimmtes Softwarepaket ad4 ist. Cross Validated betrachtet keine Fragen oder Antworten, bei denen es ausschließlich um Code geht. In Stack Overflow werden Softwareprobleme behoben.

Michael R. Chernick