Jak obliczyć główne składniki obrócone varimax w R?

Uruchomiłem PCA na 25 zmiennych i wybrałem 7 najlepszych komputerów za pomocą prcomp.

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

Następnie wykonałem obrót varimax na tych elementach.

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

A teraz chcę varimax obrócić dane obrócone PCA (ponieważ nie jest to część obiektu varimax - tylko macierz obciążeń i macierz obrotu). Przeczytałem, że w tym celu mnożymy transpozycję macierzy rotacji przez transpozycję danych, więc zrobiłbym to:

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

Ale to nie ma sensu, ponieważ powyższe wymiary macierzy transponują odpowiednio i więc pozostanie mi macierz składająca się tylko z wierszy, a nie wierszy ... czy ktoś wie co robię źle tutaj lub jaka powinna być moja ostatnia linia? Czy po prostu muszę później dokonać transpozycji?7 × 16933 7 $7\times 7$$7 \times 16933$$7$$16933$

„Obroty” to podejście opracowane w analizie czynnikowej; tam obroty (takie jak np. varimax) są stosowane do ładunków , a nie do wektorów własnych macierzy kowariancji. Ładunki są wektorami własnymi skalowanymi przez pierwiastki kwadratowe odpowiednich wartości własnych. Po rotacji varimax wektory ładujące nie są już ortogonalne (nawet jeśli obrót nazywany jest „ortogonalnym”), więc nie można po prostu obliczyć rzutów ortogonalnych danych na obrócone kierunki ładowania.

Odpowiedź FTusella zakłada, że ​​rotacja varimax jest stosowana do wektorów własnych (nie do obciążeń). To byłoby dość niekonwencjonalne. Szczegółowe informacje znajdują się na moim szczegółowym koncie PCA + varimax: Czy po PCA następuje obrót (np. Varimax) nadal PCA? W skrócie, jeśli spojrzymy na SVD macierzy danych , to obrócenie ładunków oznacza wstawienie dla niektórych macierzy rotacji w następujący sposób: R R ⊤ R X = ( U R ) ( R ⊤ S V ⊤ ) $X=USV^\top$$RR^\top$$R$$X=(UR)(R^\top SV^\top).$

Jeśli obrót zostanie zastosowany do obciążeń (jak to zwykle bywa), istnieją co najmniej trzy proste sposoby obliczania komputerów z rotacją varimax w R:

1. Są one łatwo dostępne za pośrednictwem funkcji psych::principal(co pokazuje, że jest to rzeczywiście standardowe podejście). Pamiętaj, że zwraca standardowe wyniki , tzn. Wszystkie komputery mają wariancję jednostkową.

2. Można ręcznie użyć varimaxfunkcji do obracania obciążeń, a następnie użyć nowych obróconych obciążeń, aby uzyskać wyniki; należy pomnożyć dane z transponowanym pseudo-odwrotnością obróconych ładunków (patrz wzory w tej odpowiedzi przez @ttnphns ). To da również znormalizowane wyniki.

3. Można użyć varimaxfunkcji, aby obrócić ładunki, a następnie użyć $rotmatmacierzy obrotu, aby obrócić znormalizowane wyniki uzyskane za pomocą prcomp.

Wszystkie trzy metody dają ten sam wynik:

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

Daje to trzy identyczne wyniki:

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

---

Uwaga:varimax funkcja w R wykorzystuje normalize = TRUE, eps = 1e-5parametry domyślnie ( patrz dokumentacja ). Można chcieć zmienić te parametry (zmniejszyć epstolerancję i zadbać o normalizację Kaiser) podczas porównywania wyników z innym oprogramowaniem, takim jak SPSS. Dziękuję @GottfriedHelms za zwrócenie mojej uwagi na to. [Uwaga: te parametry działają po przekazaniu do varimaxfunkcji, ale nie działają po przekazaniu do psych::principalfunkcji. Wygląda to na błąd, który zostanie naprawiony.]

Musisz użyć matrycy $loadings, a nie $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

Macierz $rotmatjest macierzą ortogonalną, która wytwarza nowe ładunki z nieobrotowych.

EDYTUJ od 12 lutego 2015 r .:

Jak słusznie wskazał poniżej @amoeba (patrz również jego / jej poprzedni post oraz inny post z @ttnphns ), ta odpowiedź jest nieprawidłowa. Rozważmy macierzy danych . Rozkład wartości w liczbie pojedynczej to gdzie ma za kolumny (znormalizowane) wektory własne . Teraz obrót jest zmianą współrzędnych i sprowadza się do zapisania powyższej równości jako: przy czym jest macierzą ortogonalną wybraną do uzyskania bliski rzadkości (maksymalny kontrast między wpisami, luźno mówiąc). Teraz, jeśli to wszystko$n\times m$$X$

$$X = USV^T$$ $V$$X'X$ $$X = (UST)(T^TV^T) = U^*V^*$$ $T$$V^*$, co nie jest, można pomnożyć powyższą równość przez aby uzyskać wyniki jako , ale oczywiście nigdy nie obracamy całego komputera. Przeciwnie, uważa się, że podzbiór , który zapewnia jeszcze godnego rank- aproksymacji , tak więc rozwiązanie jest obrócony gdzie teraz jest macierzą . Nie możemy już po prostu pomnożyć przez transpozycję , ale raczej musimy skorzystać z jednego z rozwiązań opisanych przez @amoeba.$V^*$$U^*$$X(V^*)^T$$k<m$$k$$X$X≈(UkSkTk)(T T k V T k )=U ∗ k V ∗ k V ∗ k k×nXV ∗ $$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^*$

Innymi słowy, zaproponowane przeze mnie rozwiązanie jest poprawne tylko w tym konkretnym przypadku, w którym byłoby bezużyteczne i bezsensowne.

Serdeczne podziękowania dla @amoeba za wyjaśnienie mi tej sprawy; Od lat żyję z tym nieporozumieniem.

Jeden punkt, w którym uwaga powyżej odbiega od użytkownika @ ameba postu jest to, że ona / on wydaje skojarzyć z w . Myślę, że w PCA bardziej powszechne jest posiadanie kolumn normy 1 i pochłanianie w wartościach głównego składnika. W rzeczywistości zwykle są one przedstawiane jako kombinacje liniowe oryginalnych (wyśrodkowanych, być może skalowanych) zmiennych podlegających . Tak czy inaczej jest do zaakceptowania, jak sądzę, i wszystko pomiędzy (jak w analizie biplot).$S$$V$$L$$V$$S$$v_i^TX$ $(i=1,\ldots,m)$$\|v_i\|=1$

DALSZA EDYCJA 12 lutego 2015 r

Jak zauważył @amoeba, mimo że jest prostokątny, zaproponowane mnie rozwiązanie może być nadal akceptowalne: dałoby macierz jednostkową i . Wszystko wydaje się więc zależeć od definicji wyników, które preferuje się. V ∗ k ( V ∗ k ) T X ( V ∗ k ) T ≈ U ∗ $V_k^*$$V_k^*(V_k^*)^T$$X(V_k^*)^T \approx U_k^*$

Szukałem rozwiązania, które działa na PCA wykonywane przy użyciu ade4 .

Poniżej znajdziesz funkcję:

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

^{Utworzono 14.01.2020 przez pakiet reprezentx (v0.3.0)}

Mam nadzieję, że to pomoże!